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1  Introduction 


Due  to  severe  resolution  requirements  resulting  in  excessive  expenditure  of  computational  resources,  direct 
numerical  simulation  (DNS)  of  turbulent  flows  is  generally  limited  to  relatively  low  Reynolds  numbers  and 
to  simple  geometric  configurations.  In  order  to  reduce  these  demands,  particularly  for  practical  applications, 
it  is  desirable  to  model  certain  aspects  of  the  turbulence  in  some  manner.  This  has  been  difficult  because 
the  large-scale  structures,  which  contain  most  of  the  turbulent  energy,  vary  considerably  from  one  flow  to 
another,  thereby  precluding  a  general  description.  In  large-eddy  simulation,  only  the  finest  structures  are 
left  under-resolved,  and  their  dominant  effect  must  be  accounted  for  through  some  other  means.  Since  fine- 
scale  structures  are  believed  to  be  homogeneous  and  possess  a  universal  character,  their  effect  may  be  more 
easily  and  reliably  modeled.  Additionally,  the  fine  structures  contain  only  a  fraction  of  the  total  turbulent 
kinetic  energy,  and  it  is  therefore  generally  assumed  that  they  may  be  accounted  for  without  unduly  affecting 
larger  turbulent  eddies.  This  may  be  done  through  use  of  an  SGS  stress  model  to  provide  the  dissipation 
supplied  by  the  fine-scale  structures,  or  by  the  dissipation  inherent  in  the  numerical  solving  procedure.  In 
the  latter  case,  care  must  be  taken  that  implicit  dissipation  of  the  computational  scheme  does  not  dominate 
that  intrinsic  to  an  explicitly  added  SGS  model  approach. 

Applications  of  LES  to  increasingly  practical  configurations  of  engineering  interest,  is  motivated  by 
the  need  to  provide  more  realistic  characterizations  of  the  complex  unsteady  and  separated  flows  that  are 
encountered  in  areas  of  aeroacoustics,  aero-optics,  fluid/structure  interactions,  and  active  flow  control.  In 
these  situations,  accurate  prediction  of  the  flowfields  requires  that  the  large-scale  dynamics  must  be  properly 
captured,  which  is  a  requirement  beyond  the  capabilities  of  traditional  Reynolds-averaged  Navier-Stokes 
(RANS)  methodology.  Particularly  for  control  applications,  the  baseline  case  that  is  sought  to  be  modified, 
may  contain  large  regions  of  highly  unsteady  separated  flow  that  are  not  in  turbulent  equilibrium.  Small- 
scale  fluid  structures  may  be  present,  and  portions  of  the  flowfield  may  also  be  transitional.  For  these 
reasons,  descriptions  by  RANS  models  are  highly  inadequate.  In  addition,  control  techniques  are  often 
predicated  upon  unsteady  forcing  which  can  generate  additional  small-scale  structures,  enhance  mixing,  and 
create  supplemental  turbulent  kinetic  energy.  Unsteady  forcing  may  also  be  used  to  perturb  unstable  shear 
layers,  or  as  a  “tripping”  mechanism  to  promote  bypass  transition.  Because  more  recently  developed  hybrid 
procedures[l,  2,  3]  that  combine  RANS  modeling  with  large-eddy  simulation  have  not  been  shown  to  be 
generally  satisfactory  for  active  flow  control  applications,  numerical  techniques  at  least  as  sophisticated  as 
LES  are  more  commonly  required  for  this  class  of  problems. 

The  purpose  of  this  paper  is  to  describe  a  computational  method  for  performing  large-eddy  simulation, 
that  has  been  developed  over  an  approximately  ten  year  period.  As  the  method  has  matured  during  that 
time,  a  number  of  active  flow  control  applications  have  been  considered,  some  of  which  are  summarized  here. 
There  are  several  features  of  the  general  approach  which  make  it  attractive  for  performing  active  flow  control 
simulations.  It  is  based  upon  an  implicit  time-marching  algorithm[4],  so  that  it  is  suitable  for  wall-bounded 
flows.  High-order  spatial  accuracy  is  achieved  by  use  of  an  implicit  compact  finite-difference  scheme [5], 
making  LES  resolution  attainable,  with  a  minimal  expenditure  of  computational  resources.  Robustness  is 
enhanced  by  employing  a  low-pass  Pad e- type  non-dispersive  spatial  filter  [6]  that  regularizes  the  solution  in 
flow  regions  where  the  computational  mesh  is  not  sufficient  to  fully  resolve  the  smallest  captured  structures. 
Due  to  the  spectral-like  dissipation  properties  of  the  filter,  it  also  serves  the  same  function  as  that  of  an 
SGS  model  without  additional  computational  expense.  In  order  to  accommodate  geometrically  complex 
configurations,  an  overset  grid  technique  is  adopted,  with  high-orcler  interpolation [7,  8]  to  maintain  spatial 
accuracy  at  overlapping  mesh  interfaces.  Details  of  these  features  are  described  in  sections  that  follow. 
To  illustrate  application  for  active  flow  control  simulations,  several  separate  computations  are  considered. 
These  consist  of  acoustic  resonance  suppression  in  supersonic  cavity  flow,  leading-edge  vortex  control  of  a 
delta  wing,  efficiency  enhancement  of  a  transitional  highly-loaded  low-pressure  turbine  blade,  and  separation 
control  of  a  wall-mounted  hump  model.  Control  techniques  represented  in  these  examples  are  comprised 
of  both  steady  and  pulsed  mass  injection  or  removal,  as  well  as  plasma-based  actuation.  Features  of  the 
flowfield  are  elucidated  for  each  case,  and  comparisons  are  made  with  baseline  situations  where  no  control 
was  enforced,  and  with  experimental  data  where  available. 
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2  The  Governing  Equations 


The  governing  fluid  equations  are  taken  as  the  unsteady  three-dimensional  compressible  unfiltered  Navier- 
Stokes  equations.  Although  these  computations  are  considered  to  be  large-eddy  simulations,  it  will  be 
subsequently  explained  why  the  unfiltered  equations  are  solved.  After  introducing  a  transformation  from 
Cartesian  coordinates  to  a  general  time-dependent  body-fitted  curvilinear  system,  the  equations  are  cast  in 
the  following  nondimensional  conservative  form 
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Here  t  is  the  time,  £,  77,  C,  the  computational  coordinates,  Q  the  vector  of  dependent  variables,  F ,  G,  H  the 
inviscid  flux  vectors,  and  Fv,  Gv,  Hv  the  viscous  flux  vectors.  A  source  vector  S  has  been  included  in  the 
formulation,  and  is  used  to  represent  the  body  force  induced  by  an  electric  field  for  examples  which  utilized 
plasma-based  control.  The  vector  of  dependent  variables  is  given  as 
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In  the  preceding  expressions,  u,v,w  are  the  Cartesian  velocity  components,  p  the  density,  p  the  pres¬ 
sure,  and  T  the  temperature.  All  length  scales  have  been  nondimensionalized  by  the  characteristic  length 
l,  and  dependent  variables  have  been  normalized  by  their  reference  values  except  for  p  which  has  been 
nondimensionalized  by  Pooii ‘^yo.  Components  of  the  stress  tensor  and  heat  flux  vector  are  expressed  as 
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The  Sutherland  law  for  the  molecular  viscosity  coefficient  /x  and  the  perfect  gas  relationship 


pT 

7  Ml 


(11) 


were  also  employed,  and  Stokes’  hypothesis  for  the  bulk  viscosity  coefficient  has  been  invoked. 

During  development  of  the  computational  procedure,  the  use  of  traditional  SGS  stress  models  was  ex¬ 
plored,  and  had  been  employed  for  one  of  the  example  problems  presented  below  (supersonic  cavity) .  In  that 
circumstance,  the  filtered  version  of  the  governing  equations  was  considered.  The  filtered  form  is  similar 
to  the  one  presented  above,  but  includes  additional  SGS  stress  and  heat  flux  terms  which  necessarily  must 
be  modeled.  Principal  effects  of  these  terms  are  embodied  in  modified  values  of  the  molecular  viscosity 
coefficient  /x  and  the  Prandtl  number  Pr,  which  are  augmented  to  include  contributions  due  to  the  SGS 
model. 
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The  Numerical  Method 


3.1  The  Time  Marching  Scheme 

Time-accurate  solutions  to  Eq.  (1)  were  obtained  numerically  by  the  implicit  approximately-factored  finite- 
difference  algorithm  of  Beam  and  Warming[4]  employing  Newton-like  subiterations, [9]  which  has  evolved  as 
an  efficient  tool  for  generating  solutions  to  a  wide  variety  of  complex  fluid  flow  problems.  The  algorithm  is 
used  to  advance  the  solution  in  time,  and  may  be  written  in  “delta”  form  as 
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In  this  expression,  Qp+1  is  the  p+1  approximation  to  Q  at  the  n+1  time  level  Qn+1,  and  A  Q  =  Qp+1—  Qp. 
For  p  =  1,  Qp  =  Qn .  Second-order-accurate  backward-implicit  time  differencing  was  used  to  obtain  temporal 
derivatives. 

The  implicit  segment  of  the  algorithm  incorporates  second-order-accurate  centered  differencing  for  all 
spatial  derivatives.  Nonlinear  artificial  dissipation  terms[10,  11]  are  also  appended  to  the  implicit  operators 
to  augment  stability,  but  for  simplicity  have  not  been  shown  explicitly  in  Eq.  (12).  Efficiency  is  enhanced 
by  solving  this  implicit  portion  of  the  factorized  equations  in  diagonalized  form.  [12]  Temporal  accuracy, 
which  can  be  degraded  by  use  of  the  diagonal  form,  is  maintained  by  utilizing  subiterations  within  a  time 
step.  This  technique  has  been  commonly  invoked  in  order  to  reduce  errors  due  to  factorization,  lineariza¬ 
tion,  diagonalization,  and  explicit  application  of  boundary  conditions.  It  is  useful  for  achieving  temporal 
accuracy  on  overset  zonal  mesh  systems,  and  for  a  domain  decomposition  implementation  on  parallel  com¬ 
puting  platforms.  Any  deterioration  of  the  solution  caused  by  use  of  artificial  dissipation  and  by  lower-order 
spatial  resolution  of  implicit  operators  is  also  reduced  by  the  procedure.  From  a  large  number  of  previous 
computations,  it  was  found  that  three  subiterations  per  time  step  were  sufficient  to  preserve  second-order 
temporal  accuracy.  Because  temporal  accuracy  is  limited  by  the  backward  approximation  in  Eq.  (12),  further 
subiterations  do  not  improve  the  solution.  For  all  of  the  examples  considered  here,  three  subiterations  per 
time  step  have  been  applied. 


3.2  The  Compact  Finite-Difference  Scheme 

The  compact  difference  scheme  employed  to  obtain  spatial  derivatives  on  the  right-hand  side  of  Eq.  (12)  is 
based  upon  the  pentadiagonal  system  of  Lele,[5]  and  is  capable  of  attaining  spectral-like  resolution.  This  is 
achieved  through  the  use  of  a  centered  implicit  difference  operator  with  a  compact  stencil,  thereby  reducing 
the  associated  discretization  error.  The  scheme  is  illustrated  here  in  one  spatial  dimension  on  a  uniformly 
spaced  mesh  for  the  general  function  <fi(x)  as 
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Equation  (13)  can  be  used  to  define  families  of  both  explicit  (a  =  (3  =  0)  and  implicit  difference  approxima¬ 
tions  by  proper  choice  of  the  coefficients  a,  (3,  a,  6,  c.  Because  all  of  these  schemes  retain  the  central-difference 
formulation,  there  is  no  dissipation  error  associated  with  any  of  them.  Following  Lele[5]  however,  it  is  useful 
to  examine  their  dispersion  characteristics.  This  is  done  by  performing  a  spatial  wave  number  analysis  on 
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the  Fourier  components  of  the  function  (f>(x).  For  the  analysis,  we  assume  4>{x)  is  periodic  on  the  interval 
0  <  x  <  L  and  let 
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and  <j)m  are  the  Fourier  coefficients.  It  is  then  convenient  to  define  the  scaled  wave  number  O  and  the  scaled 
coordinate  S  by 

O  =  2nmAx  and  s  =  respectively,  (15) 


so 
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and  to  denote  d<f>/dS  =  (j).  From  Eq.  (16),  it  follows  that 
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and  upon  comparing  Eq.  (16)  with  Eq.  (17)  it  is  seen  that  the  exact  derivative  of  <f>  generates  Fourier 
coefficients  that  are  simply  related  to  those  of  the  original  function  by  the  expression 
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By  direct  substitution  of  Eq.  (16)  into  the  difference  formula  Eq.  (13),  the  numerical  approximation  to  the 
derivative  (</>)n  results  in  Fourier  coefficients  with  the  relationship 


(0m) n  =  0m. 

where  fi'  is  a  modified  wave  number  satisfying  the  expression 

asin(fl)  +  (6/2)  sin(2fl)  +  (c/3)  sin(3fl) 


Q'  = 
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Shown  in  Fig.  1  is  a  plot  of  the  modified  wave  number  Q'  as  a  function  of  the  wave  number  f 1  for 
a  number  of  explicit  and  implicit  schemes.  For  the  function  <f>,  2tt/Q  indicates  the  number  of  points  per 
period  that  are  represented  locally.  The  exact  solution  is  given  by  Cl'  =  O.  It  is  seen  that  the  compact 
implicit  schemes  are  able  to  better  approximate  the  exact  solution  at  much  higher  wave  numbers.  This  is 
equivalent  to  requiring  fewer  points  per  period  in  order  to  achieve  the  same  resolution.  Lele[5]  defined  a 
resolving  efficiency  for  comparing  the  various  methods.  From  this  criterion,  it  was  shown  that  the  sixth-order 
tridiagonal  (/?  =  c  =  0)  subset  of  Eq.  (13)  has  5-10  times  better  resolution  than  the  traditional  2nd-order 
explicit  approach  (a  =  /3  =  6  =  c=0,  a  =  1.0).  It  was  further  shown,  as  can  be  inferred  from  Fig.  1,  that 
the  resolving  power  of  the  lOth-orcler  pentadiagonal  formulation  was  not  appreciably  better  than  that  of 
the  sixth-order  tridiagonal  subset.  We  therefore  restrict  our  spatial  differencing  scheme  to  a  4th-order  or 
6th-order  tridiagonal  implicit  approximation,  where 


a  =  1/4,  (3  =  0,  a  =  3/2,  b  =  c  =  0  for  the  4th-order  scheme 
and  a  =1/3,  (3  =  0,  a  =  14/9,  6=1/9,  c=0  for  the  6th-order  scheme. 


(21) 

(22) 


Based  upon  our  experience,  the  tridiagonal  subset  of  Eq.  (13)  increases  the  computational  time  by  about 
a  factor  of  two  over  that  of  a  standard  2nd-order  explicit  scheme  for  solution  of  Eq.  (1).  But  because  of 
superior  resolving  capability,  fewer  computational  resources  need  be  expended  with  the  high-orcler  method, 
than  are  required  with  the  standard  approach,  in  order  to  attain  the  same  level  of  resolution.  Solution 
of  the  tridiagonal  system  of  Eq.  (13)  is  about  50%  computationally  less  expensive  than  the  pentadiagonal 
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counterpart.  Thus  the  4th-order  and  6th-order  compact  difference  schemes  provide  a  somewhat  optimal 
balance  between  efficiency  and  accuracy. 

The  compact  difference  schemes  described  above  are  used  to  obtain  spatial  derivatives  of  any  scalar 
variable,  such  as  a  flux  component,  flow  quantity,  or  metric  coefficient  appearing  on  the  right-hand  side  of 
Eq.  (12).  Derivatives  of  inviscid  fluxes  are  computed  by  forming  flux  quantities  at  grid  points,  and  then 
applying  the  above  formulas  to  each  component.  Viscous  derivatives  are  obtained  by  first  computing  high- 
order  derivatives  of  primitive  variables.  Components  of  the  viscous  fluxes  are  then  constructed,  and  the 
compact  difference  scheme  is  applied  a  second  time.  Although  this  technique  is  not  formally  as  accurate  as 
a  high-order  scheme  employed  directly  for  evaluation  of  second  derivatives,  it  requires  less  computational 
effort.  It  has  also  been  demonstrated  to  produce  accurate  and  stable  results  for  for  LES  computations.  [13] 
A  description  of  high-order  one-sided  difference  formulas  for  use  near  boundaries  may  be  found  in  Ref.  [14]. 

The  compact  difference  scheme  is  also  used  to  evaluate  metric  coefficients  and  the  Jacobian  of  the  coor¬ 
dinate  transformation  indicated  in  Eq.  (12).  This  is  done  for  example,  in  the  relationship 

ix  =  yvzc  -  y^zy  (23) 

by  applying  difference  formulas  to  the  analytic  equivalent  conservative  form 

ix  =  {yvz)c  -  (ycz)r,.  (24) 

Use  of  Eq.  (24)  follows  the  development  from  Refs.  [15]  and  [16]  which  are  based  upon  the  treatment  by 
Thomas  and  Lombard[17]  for  low-order  methods.  The  technique  is  used  to  preserve  the  metric  identities 

+  (Vx)r)  +  (Cx)c  =  (£y)j  +  ( rly)v  +  (Cy)c  =  (£z){  +  ('lh)v  +  (Cz)c  =  0  (25) 

which  were  implicitly  invoked  for  the  derivation  of  Eq.  (12).  Reference  [18]  gives  additional  details  about 
the  treatment  of  time-dependent  metric  coefficients  and  the  Jacobian  for  deforming  mesh  applications. 

3.3  The  Compact  Filtering  Scheme 

As  noted  previously,  the  compact  central  difference  approximation  is  nondissipative,  and  therefore  can  be 
susceptible  to  numerical  instabilities  which  arise  due  to  unrestricted  growth  of  high-frequency  spatial  modes. 
Such  instabilities  originate  from  several  sources,  that  include  mesh  nonuniformities,  approximate  boundary 
conditions,  and  nonlinear  flow  behavior.  In  order  to  extend  the  compact  discretization  to  practical  appli¬ 
cations,  a  high-order  low-pass  Pade-type  non-dispersive  spatial  filtering  technique[13,  6]  is  incorporated  as 
part  of  the  numerical  methodology.  This  low-pass  filter  provides  dissipation  at  high  spatial  wave  numbers, 
only  where  the  resolution  already  exhibits  significant  dispersion  error. 

A  general  expression  for  a  family  of  implicit  filters  applied  to  the  function  (j>  may  be  written  as 

N 

Pf4>i- 2  +  1  +  4>i  +  af(t>i+ 1  +  Pf4>i+ 2  =  ~Jr  {4>i-n  +  4>i+n)  (26) 

n—0 

where  <f>  is  the  filtered  value  of  <f>.  Equation  (26)  can  be  used  to  define  unique  filter  representations  with 
accuracy  of  order  27V+6.  These  filters  provide  no  dispersion,  and  once  more  following  Lele[5],  their  dissipation 
properties  may  characterized  from  a  wave  number  analysis.  Assuming  again  that  4>  is  periodic  and  may  be 
expressed  in  terms  of  the  Fourier  series  given  by  Eq.  (16),  it  follows  from  Eq.  (26)  that 


m=+M/ 2 

4>  =  T(U)0m  exp(fUS') 

m=-M/2 


where  T(f2) 
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Equation  (27)  illustrates  that  Fourier  components  of  the  filtered  function  <f>m  are  related  to  those  of  the 
original  function  by 

(29) 

where  T(f2)  is  the  spectral  transfer  function.  We  find  it  practical  to  restrict  consideration  to  tridiagonal 
subsets  of  Eq.(26)  (/?/  =  0).  Although  it  is  then  possible  to  achieve  a  unique  description  of  order  2 TV  +  4 
using  a  2 N  +  1  stencil,  we  choose  to  limit  the  order  to  2 TV.  This  allows  freedom  to  enforce  two  additional 
constraints  upon  the  filter  operator  through  the  choice  of  the  coefficients  a/,a„’ s.  One  of  these  is  supplied 
by  requiring 

T(tt)  =  0  (30) 

which  fully  damps  all  contributions  with  resolution  less  than  two  points  per  period.  The  other  condition 
allows  the  implicit  coefficient  a/  to  remain  a  free  parameter,  which  may  be  adjusted  for  specific  applications. 
For  proper  behavior  of  the  transfer  function  T(f l),  the  adjustable  coefficient  a/  must  lie  in  the  range  —0.5  < 
ctf  <  0.5,  where  higher  values  correspond  to  less  dissipative  filters.  Explicit  filter  formulas  are  obtained  with 
cif  =  0.  On  uniform  meshes,  these  symmetric  filters  are  non-dispersive  (T(f2)  is  real),  do  not  amplify  waves 
(T(fi)  <  1),  preserve  constant  functions  (T(f2)  =  1),  and  preclude  odd-even  mode  decoupling  (T(7t)  =  0). 

Because  the  finite-difference  approximation  is  limited  to  at  most  sixth-order  accuracy,  we  consider  filters 
no  higher  than  lOth-order,  having  an  11  point  stencil,  so  that  our  filtering  formula  is  given  by 

N 

1  +  <j>i  +  Ctf(j>i+1  =  ( 4>i-n  +  <j>i+n)  N  <5.  (31) 

n— 0 

It  is  required  that  accuracy  of  the  filter  should  be  at  least  two  orders  higher  than  that  of  the  differencing 
scheme  employed  for  any  large-eddy  simulation.  Presented  in  Fig.  2  is  a  plot  of  the  transfer  function  T (fi) 
for  several  different  filters.  The  exact  result  corresponds  to  the  unfiltered  form  with  T(f 1)  =  1.  As  expected, 
lower-order  filters  provide  more  dissipation  at  lower  wave  numbers.  Also,  for  filters  of  the  same  order,  explicit 
formulations  are  more  dissipative  than  implicit  ones  (see  6th-order  results  in  Fig.  2).  The  effect  of  varying 
the  implicit  coefficient  a/  for  the  lOth-order  tridiagonal  subset  is  illustrated  in  Fig.  3.  With  a/  =  0.49,  the 
spectral-like  behavior  is  demonstrated. 

Unlike  numerical  methods  which  employ  the  use  of  explicitly  added  artificial  dissipation,  thereby  mod¬ 
ifying  the  original  governing  equations,  the  filtering  operation  is  a  post  processing  technique.  It  is  applied 
to  the  evolving  solution  in  order  to  regularize  features  that  are  captured  but  poorly  resolved.  Filtering  is 
carried  out  on  the  conserved  variables  successively  along  each  of  the  transformed  coordinate  directions.  For 
the  example  applications  subsequently  to  be  described,  filtering  was  imposed  once  following  every  subitera¬ 
tion.  In  general  however,  it  may  be  invoked  repeatedly  or  less  frequently  depending  on  the  nature  of  specific 
problems.  Numerical  values  for  the  explicit  coefficients  an’s  as  a  function  of  the  implicit  coefficient  a /  may 
be  found  in  Ref.  [14] . 

A  tenth-order  central  filter  having  an  11-point  stencil  cannot  be  applied  within  a  distance  of  five  grid 
points  from  a  computational  boundary.  In  this  region,  a  centered  stencil  may  be  retained  only  if  the  orcler-of- 
accuracy  of  the  filter  is  reduced,  and/or  one-sided  filter  formulations  are  employed.  As  opposed  to  the  central 
formulation,  one-sided  filters  may  generate  amplification  in  certain  ranges  of  wave  number.  Amplification 
levels  become  larger  as  the  order  of  the  filter  is  increased,  corresponding  to  a  more  one-sided  filter.  The 
region  over  which  this  amplification  occurs  however,  shifts  to  higher  wave  numbers  with  increasing  order, 
thus  improving  spectral  behavior  in  the  resolved  wave  space  region.  Spurious  waves  that  may  be  amplified 
by  high-order  one-sided  boundary  filters  tend  to  be  removed  by  interior  filters  as  they  propagate  away 
from  boundaries.  Increasing  the  value  of  the  filter  coefficient  a /  for  near-boundary  points  can  also  reduce 
amplification  associated  with  one-sided  filters.  Filtering  strategies  for  the  treatment  of  near-boundary  points 
appear  in  Refs.  [13]  and  [15].  The  impact  of  filtering  on  the  accuracy  and  stability  of  the  high-order  approach 
has  been  investigated  by  Visbal  and  Gaitonde[13,  15,  16]  for  a  number  of  flow  situations,  which  address 
nonuniform  grids,  approximate  boundary  treatments,  and  nonlinear  governing  equations. 

3.4  The  Overset  Grid  Approach 

An  overset  grid  approach  [19]  is  employed  in  order  to  provide  flexibility  for  treating  complex  geometric 
configurations.  This  allows  structured  grid  formulations  of  the  compact  differencing  and  filtering  schemes 
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to  be  adapted  for  practical  simulations.  In  addition,  the  technique  also  allows  local  grid  refinement  which 
may  be  necessary  to  capture  physical  details  arising  from  a  large-eddy  simulation.  It  is  particularly  useful 
for  consideration  of  active  flow  control,  that  may  require  modeling  of  actuation  devices,  or  may  locally 
generate  small-scale  fluid  structures.  Furthermore,  local  refinement  can  be  used  to  reduce  overall  resource 
requirements  for  computationally  demanding  simulations. 

Although  the  overset  grid  approach  is  fundamentally  motivated  by  its  ability  to  describe  geometric  com¬ 
plexity,  it  is  also  employed  to  enable  domain  decomposition  for  processing  on  parallel  computing  platforms. 
Individual  unique  mesh  entities  may  be  sub-divided  into  any  number  of  smaller  grid  systems,  each  of  which 
are  then  operated  on  by  a  single  processor.  To  maintain  high-order  accuracy,  all  grid  systems  are  required  to 
overlap  with  their  adjoining  neighbors.  This  is  true  for  both  sub-divided  meshes  and  geometrically  distinct 
grids.  It  has  been  shown  in  previous  studies[20],  that  an  overlap  of  five-planes  in  the  region  between  respec¬ 
tive  meshes,  is  sufficient  to  maintain  the  interior  high-order  differencing  and  filtering  schemes.  The  overlap 
region  consists  of  two  levels  (contiguous  planes)  of  donor/receiver  planes  of  grid  points  which  are  overset  into 
adjacent  domains  on  each  side  of  the  overlap.  This  arrangement  is  typically  employed  for  decomposed  grids 
with  coincident  mesh  points,  but  a  somewhat  larger  overlap  may  necessarily  be  utilized  for  general  overset 
domains  due  to  variations  in  grid  topology. 

An  example  of  the  five-point  overlap  for  domain  decomposition  is  shown  in  Fig.  4.  The  figure  illustrates 
details  of  the  overlap  for  the  simple  schematic  representation  of  a  vortex  convecting  between  two  sub-domains 
separated  by  a  vertical  interface.  The  vortex  located  in  mesh  1,  is  traveling  from  left  to  right.  It  passes 
through  the  five  point  overlap  and  into  mesh  2.  Information  is  transferred  between  the  two  sub-domains, 
which  are  distributed  on  different  processors,  via  the  five-point  overlap  using  message-passing  interface 
(MPI)  communications.  [21]  The  expanded  view  of  the  overlap  region  in  the  figure  shows  the  two  sets  of  five 
vertical  lines,  each  denoted  by  its  I  index,  which  are  used  for  communication  between  the  grids.  Although 
the  overlap  points  are  coincident,  they  have  been  shown  slightly  staggered  for  clarity.  Data  is  exchanged 
between  adjacent  subdomains  at  the  end  of  each  sub-iteration  of  the  time-marching  scheme.  The  values  of 
flow  variables  at  points  1  and  2  of  mesh  2  are  set  to  be  equal  to  those  of  the  corresponding  updated  values 
at  points  IL  —  4  and  IL  —  3  of  mesh  1.  Similarly,  reciprocal  information  is  transferred  through  points  4  and 
5  of  mesh  2  which  “inject”  the  solution  into  points  IL  —  1  and  IL  of  mesh  1.  Point  3  in  mesh  1  and  IL  —  2 
in  mesh  2  are  solved  independently,  and  not  updated  with  the  solution  from  the  adjacent  domain.  Use  of 
these  dual  solutions,  facilitates  detecting  any  disparity  between  results  in  the  respective  domains. 

Connectivity  between  donor/receiver  pairs  is  establish  as  a  pre-processing  step  prior  to  computation,  using 
automated  software.  For  each  receiver  grid  point,  the  PEGASUS [22]  utility  is  used  to  identify  a  donating 
stencil.  PEGASUS  also  defines  second-order  accurate  interpolation  coefficients.  In  the  case  of  sub-divided 
grids  with  coincident  points  in  the  overlap  region,  these  coefficients  degenerate  to  “direct  injection”  so  that 
there  is  no  interpolation  error.  When  the  overlap  region  consists  of  non-coincident  grid  points,  high-orcler 
interpolation  must  be  implemented  in  order  to  maintain  the  overall  accuracy  of  the  numerical  methodology. 
This  is  achieved  with  a  second  pre-processing  utility  BELLERO[23],  that  is  employed  to  modify  second-order 
interpolation  stencils  in  order  to  accommodate  a  more  accurate  approximation.  If  (j>v  denotes  the  interpolated 
value  of  the  flowfield  variable  <f>  at  grid  point  p  using  known  information  of  <f>  on  another  grid  at  the  donor 
stencil  defined  by  point  ( Ip ,  Jp,  Kp),  then  the  high-order  explicit  Lagrangian  formula 

a—  1  a—  1  <7  —  1 

&p  =  X!  ^2  X!  RitfRkfap+i.Jp+LKp+k  (32) 

2=0  j= 0  k—0 

is  applied.  In  Eq.(32),  a  is  the  formal  order-of-accuracy  of  the  interpolated  approximation  and  R] .  R‘j .  Rf 
are  coefficients  given  by  the  analytic  expressions 
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where  0  <  i,j,k  <  a  —  l  and  Ai,  A2,A3  are  interpolation  offsets  in  the  directions,  representing 

distances  from  the  base  donor  point  (Ip,  Jp,  Kp)  to  the  interpolation  point  in  the  computational  space  of  the 
donor  grid  (0  <  Ai,A2,A3  <  a  —  1).  The  base  donor  point  location  (. Ip,Jp,Kp )  and  interpolation  offsets 
(Ai,  A2,  A3)  are  obtained  from  PEGASUS/BELLERO. 

Once  donating  stencils  and  interpolation  coefficients  are  defined,  inter-node  communication  among  various 
processors  is  established  through  MPI  library  routines,  [21]  which  are  used  to  transfer  information  between 
the  various  grids  at  domain  boundaries,  as  previously  mentioned.  We  note  that  the  BELLERO  utility  has  a 
number  of  additional  useful  capabilities.  It  can  be  employed  as  “stand  alone”  software  (without  PEGASUS) 
for  the  partitioning  of  single  grid  systems,  where  domain  decomposition  is  enforced  only  for  parallel  processing 
and  no  true  interpolation  is  required.  BELLERO  can  also  generate  connectivity  for  application  of  periodic 
boundary  conditions,  and  may  be  employed  for  the  treatment  of  “holes”  or  “blanked  out”  regions.  The 
flowfield  in  these  “holes”  or  “blanked  out”  regions  are  then  described  by  overset  background  grids,  which  is 
a  commonly  employed  technique  in  overset  methods. 


3.5  The  LES  Approach 

With  a  traditional  LES  approach,  physical  dissipation  at  the  Kolmogorov  scale  is  not  represented.  For 
spatially  nondissipative  numerical  schemes,  without  use  of  SGS  models,  this  leads  to  an  accumulation  of 
energy  at  high  mesh  wave  numbers,  and  ultimately  to  numerical  instability.  Explicitly  added  SGS  models 
are  then  employed  as  a  means  to  dissipate  this  energy.  In  the  present  methodology,  the  effect  of  the  smallest 
fluid  structures  is  accounted  for  by  an  implicit  large-eddy  simulation  (ILES)  technique,  which  has  been 
successfully  utilized  for  a  number  of  turbulent  and  transitional  computations,  some  of  which  will  subsequently 
be  described.  The  ILES  approach  was  first  introduced  by  Visbal  et  al.  [18,  24]  as  a  formal  alternative 
to  conventional  methodologies,  and  is  predicated  upon  the  high-order  compact  differencing  and  low-pass 
spatial  filtering  schemes,  without  the  inclusion  of  additional  SGS  modeling.  This  technique  is  similar  to 
monotonically  integrated  large-eddy  simulation  (MILES)  [25]  in  that  it  relies  upon  the  numerical  solving 
procedure  to  provide  the  dissipation  that  is  typically  supplied  by  traditional  SGS  models.  Unlike  MILES 
however,  dissipation  is  contributed  only  at  high  spatial  wavenumbers  where  the  solution  is  poorly  resolved, 
by  the  aforementioned  high-order  Pade-type  low-pass  filter.  This  allows  a  mechanism  for  the  turbulence 
energy  to  be  dissipated  at  scales  that  cannot  be  accurately  resolved  on  a  given  mesh  system,  in  a  fashion 
similar  to  subgrid  modeling.  For  purely  laminar  flows,  filtering  may  be  required  to  maintain  numerical 
stability  and  preclude  a  transfer  of  energy  to  high-frequency  spatial  modes  due  to  spurious  numerical  events. 
The  ILES  methodology  thereby  permits  a  seamless  transition  from  large-eddy  simulation  to  direct  numerical 
simulation  as  the  resolution  is  increased.  In  the  ILES  approach,  the  unfiltered  governing  equations  may  be 
employed,  and  the  computational  expense  of  evaluating  subgrid  models,  which  can  be  substantial,  is  avoided. 
This  procedure  also  enables  the  unified  simulation  of  flowfields  where  laminar,  transitional,  and  turbulent 
regions  simultaneously  coexist. 

It  should  also  be  noted  that  the  ILES  technique  may  be  interpreted  as  an  approximate  deconvolution 
SGS  model[26],  which  is  based  upon  a  truncated  series  expansion  of  the  inverse  filter  operator  for  the  unfil¬ 
tered  flowfield  equations.  Mathew  et  al.  [27]  have  shown  that  filtering  provides  a  mathematically  consistent 
approximation  of  unresolved  terms  arising  from  any  type  of  nonlinearity.  Filtering  regularizes  the  solution, 
and  generates  virtual  subgrid  model  terms  that  are  equivalent  to  those  of  approximate  deconvolution. 


3.6  Validation  of  the  Numerical  Method 

The  aforementioned  features  of  the  numerical  algorithm  are  embodied  in  a  parallel  version  of  the  time- 
accurate  three-dimensional  computer  code  FDL3DI,[14]  which  has  proven  to  be  reliable  for  steady  and 
unsteady  fluid  flow  problems,  including  vortex  breakdown, [28,  29]  transitional  wall  jets, [30]  synthetic  jet 
actuators,  [31]  roughness  elements, [32]  plasma  flows, [33,  34,  35,  36]  and  direct  numerical  and  large-eddy 
simulations  of  subsonic [37,  38]  and  supersonic  flowfields.  [39,  40]  The  procedure  employed  for  parallelization 
has  been  demonstrated  to  be  portable,  due  to  the  use  of  the  standardized  MPI  library.  FDL3DI  has  been 
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utilized  successfully  on  a  number  of  different  computing  platforms,  and  has  maintained  a  parallel  efficiency 
of  90%  for  up  to  320  processors. 

During  development  of  this  code,  many  aspects  of  the  above  cited  numerical  procedure  were  largely 
validated  for  a  number  of  canonical  and  fundamental  fluid  flow  problems.  Among  these  are  the  simulation 
of  decaying  compressible  isotropic  turbulence[18,  38,  41],  for  which  it  was  shown  that  better  results  were 
realized  with  use  of  the  compact  differencing/filtering  scheme  than  could  be  obtained  with  either  a  constant 
coefficient [42]  or  dynamic  Smagorinsky[43]  SGS  model[18].  For  turbulent  channel  flows[18,  38,  8,  41],  the  use 
of  local  grid  refinement  indicated  that  high  accuracy  was  achieved  with  a  minimal  number  of  grid  points,  if 
the  mesh  system  was  distributed  properly.  It  was  also  found[41]  that  the  numerical  approach  attained  DNS 
level  resolution  for  high- Reynolds  channels. 

Simulations  for  subsonic  transitional  flow  past  a  circular  cylinder [38,  8,  41]  again  illustrated  the  benefit 
of  local  grid  refinement,  and  properly  captured  the  correct  experimentally  measured  behavior,  which  was 
not  possible  when  SGS  models  were  employed.  Computations  for  a  spatially  evolving  supersonic  flat-plate 
boundary  layer  [44]  also  compared  well  with  experimental  data,  and  performed  better  than  simulations  that 
utilized  SGS  models. 
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4  Empirical  Plasma  Model 

Two  of  the  computational  examples  to  be  considered  subsequently  (delta  wing,  low-pressure  turbine  blade), 
utilized  plasma-based  actuation  to  administer  active  flow  control.  The  following  section  describes  the  ap¬ 
proach  for  simulating  this  situation.  Many  quantitative  aspects  of  the  fundamental  processes  governing 
plasma/fluid  interactions  remain  unknown  or  computationally  prohibitive,  particularly  for  transitional  and 
turbulent  flows.  These  circumstances  have  given  rise  to  the  development  of  a  wide  spectrum  of  models  with 
varying  degrees  of  sophistication,  that  may  be  employed  for  more  practical  computations.  Among  the  sim¬ 
plified  methods  focused  specifically  on  discharge/fluid  coupling  is  that  of  Roth  et  al.,  [45,  46]  who  associated 
transfer  of  momentum  from  ions  to  neutral  particles  based  upon  the  gradient  of  electric  pressure.  A  more 
refined  approach,  suitable  for  coupling  with  fluid  response  was  an  empirical  model  proposed  by  Shyy  et  al., 
[47]  using  separate  estimates  for  the  charge  distribution  and  electric  field.  Known  plasma  physics  parameters 
were  linked  to  experimental  data.  This  representation  has  been  successfully  employed  for  several  previous 
simulations  of  plasma-controlled  flows,  [33,  34,  35,  36]  and  was  also  adopted  in  the  present  examples. 

A  schematic  representation  of  a  typical  single  asymmetric  dielectric-barrier-discharge  (DBD)  plasma  ac¬ 
tuator  is  depicted  in  Fig.  5.  The  actuator  consists  of  two  electrodes  that  are  separated  by  a  thin  dielectric 
insulator,  and  mounted  on  a  body  surface.  An  oscillating  voltage,  in  the  10-15  kHz  frequency  range,  is 
applied  to  the  electrodes,  developing  an  electric  field  about  the  actuator.  When  the  imposed  voltage  is 
sufficiently  high,  the  dielectric  produces  a  barrier  discharge,  that  weakly  ionizes  the  surrounding  gas.  Mo¬ 
mentum  acquired  by  the  resulting  charged  particles  from  the  electric  field,  is  transferred  to  the  primary 
neutral  molecules  by  a  combination  of  electrodynamic  body  forces  and  poorly  understood  complex  colli- 
sional  interactions.  Because  the  bulk  fluid  cannot  respond  rapidly  to  the  high  frequency  alternating  voltage, 
the  dominant  effect  of  actuation  is  to  impose  a  time-mean  electric  field  on  the  external  flow.  In  the  numerical 
simulation  of  plasma-based  control,  the  entire  process  was  modeled  as  a  body  force  vector  acting  on  the  net 
fluid  adjacent  to  the  actuator,  which  produces  a  flow  velocity. 

The  model  for  the  geometric  extent  of  the  plasma  field  generated  by  such  an  actuator  is  indicated  in 
Fig.  6.  The  triangular  region  defined  by  the  line  segments  AB ,  BC ,  and  AC  constitutes  the  plasma  boundary. 
Outside  of  this  region  the  electric  field  is  not  considered  strong  enough  to  ionize  the  air.  [47]  The  electric  field 
has  its  maximum  value  on  segment  AC,  and  varies  linear  within  ABC.  The  peak  value  of  the  electric  field 
is  obtained  from  the  applied  voltage  and  the  spacing  between  the  electrodes.  Along  the  segment  AB,  the 
electric  field  diminishes  to  its  threshold  value,  which  was  taken  as  30  kV/cm.  [47]  The  electric  body  force  is 
equal  to  qcE  and  provides  coupling  from  the  plasma  to  the  fluid,  resulting  in  the  source  vector  S  appearing  in 
Eq.  (1).  The  direction  of  the  force  vector  depends  upon  the  ratio  AC  I  BC,  but  for  the  simulations  presented 
here,  the  vector  was  assumed  to  be  tangential  to  the  actuator  surface.  Within  the  region  ABC,  the  charge 
density  qc  is  taken  to  be  constant.  The  plasma  scale  parameter  Dc  arises  from  nondimensionalization  of  the 
governing  equations,  and  represents  the  ratio  of  the  electrical  force  of  the  plasma  to  the  inertial  force  of  the 
fluid. 

Some  specific  details  of  the  plasma  model  incorporated  in  the  example  simulations,  were  specified  sim¬ 
ilar  to  those  of  the  original  experiment  of  Shyy  et  al.[47]  Referring  to  Fig.  6,  the  distances  BC  and  AC 
nondimensionalized  by  the  characteristic  length  l  (delta  wing  root  chord,  turbine  blade  chord)  were  taken  as 
BC  —  0.018  and  AC  =  0.024  for  the  delta  wing  and  BC  =  0.0125  and  AC  =  0.0250  for  the  turbine  blade 
respectively.  For  the  purposes  of  these  computations,  it  was  assumed  that  actuators  were  mounted  flush 
with  the  configuration  surface  and  did  not  protrude  above  it.  Due  to  empiricism  of  the  formulation,  there  is 
some  ambiguity  regarding  the  value  of  the  scale  parameter  Dc. 

Operationally,  DBD  actuators  are  inherently  unsteady  devices.  Within  the  context  of  the  empirical 
model  however,  the  body  force  imposed  on  the  fluid  is  assumed  to  be  steady,  owing  to  the  high  frequency 
of  the  applied  voltage  (typically  10-15  kHz).  These  devices  may  also  be  operated  in  a  pulsed  manner  as 
described  by  Corke  and  Post,  [48]  thereby  reducing  total  power  consumption.  The  pulsed  mode  of  operation 
also  introduces  low-frequency  forcing  to  the  flow,  which  may  be  more  receptive  to  control,  and  offers  the 
potential  of  improved  effectiveness.  These  actuators  have  no  moving  parts,  can  be  surface  conforming,  and 
provide  on  demand  control  with  low  power  utilization. 
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Table  1:  Cavity  flow  conditions 


Re 

Reso 

Re-go 

present  LES 

2.00  x  105 

6,774 

667 

experiment [49,  50,  51,  52] 

3.01  x  106 

51,700 

5,030 

5  Acoustic  Suppression  of  Supersonic  Cavity  Flow 

High  speed  flows  over  open  cavities  produce  complex  unsteady  interactions,  which  are  characterized  by 
a  severe  acoustic  environment.  At  high  Reynolds  numbers,  such  flowfields  are  comprised  of  both  broad¬ 
band  small-scale  fluctuations  typical  of  turbulent  shear  layers,  and  discrete  resonance  whose  frequency  and 
amplitude  depend  upon  the  cavity  geometry  and  external  flow  conditions.  While  these  phenomena  are  of 
fundamental  physical  interest,  they  also  represent  a  number  of  significant  concerns  for  aerospace  applications. 
In  the  practical  situation  of  an  aircraft  weapons  bay,  aerodynamic  performance  or  stability  may  be  adversely 
affected,  structural  loading  may  become  excessive,  and  sensitive  instrumentation  may  be  damaged.  Acoustic 
resonance  can  also  pose  a  threat  to  the  safe  release  and  accurate  delivery  of  weapons  systems  stored  within 
the  cavity. 

Because  of  the  intricate  nature  and  practical  significance  of  supersonic  cavity  flows,  numerous  experi¬ 
mental  and  numerical  investigations  have  been  conducted  in  order  to  understand  their  underlying  physical 
behavior.  A  summary  of  these  may  be  found  in  Ref.  [40].  For  the  simulations  presented  here,  we  consider 
a  deep  cavity  where  the  length  to  depth  ratio  is  less  than  9.  In  this  situation,  the  boundary  layer  ahead 
of  the  upstream  lip  forms  into  a  free  shear  layer,  which  then  flows  over  the  cavity  and  impinges  upon  the 
rear  bulkhead.  The  undulating  shear  layer  generates  strong  compression  waves  both  external  and  internal  to 
the  cavity,  and  results  in  periodic  addition  or  removal  of  mass  from  the  cavity  at  the  downstream  bulkhead, 
thereby  producing  a  self-sustained  fluid  oscillation. 

The  purpose  of  these  simulations  was  to  numerically  reproduce  the  turbulent  supersonic  flow  past  an 
open  rectangular  cavity.  In  addition  to  the  baseline  flow  without  control,  high-frequency  forcing,  via  pulsed 
mass  injection  upstream  of  the  forward  cavity  lip,  was  implemented  numerically  in  order  to  investigate 
the  ability  of  large-eddy  simulation  to  predict  acoustic  resonance  suppression.  For  both  unsuppressed  and 
suppressed  cases,  experimental  data  was  available  for  comparative  purposes  in  the  form  of  instantaneous 
pressure  measurements[49,  50,  51,  52]  on  interior  cavity  surfaces. 

The  cavity  geometric  configuration  is  represented  schematically  in  Fig.  7,  where  the  length  to  depth 
ratio  l/d  =  5.0  and  the  freestream  Mach  number  is  1.19.  The  cavity  length  l  was  used  as  the  reference 
distance  for  nondimensionalization.  Conditions  at  which  experimental  measurements  were  taken[49,  50, 
51,  52]  correspond  to  a  Reynolds  number  of  3.01  x  106  based  upon  the  cavity  length  l.  Because  it  would 
not  be  possible  to  numerically  resolve  fine-scale  turbulent  structures  at  the  experimental  Reynolds  number, 
simulations  were  carried  out  for  Re  =  2.0  x  105.  Reference  quantities  in  terms  of  the  incoming  boundary-layer 
parameters  are  provided  in  Table  1. 

As  previously  noted,  these  large-eddy  simulations  were  carried  out  using  the  dynamic  Smagorinsky  SGS 
model[43].  While  an  entire  description  of  the  equations  governing  this  model  is  not  crucial  to  this  example, 
complete  details  are  documented  in  Ref.  [53] .  The  numerical  scheme  for  the  computations  utilized  a  fourth- 
order  compact  finite-difference  approximation  and  a  sixth-order  spatial  filter.  A  mesh  system  consisting  of 
two  discrete  domains  was  employed  to  represent  the  regions  exterior  and  interior  to  the  cavity  flowfield.  These 
consisted  of  (725  x  225  x  101)  and  (351  x  121  x  101)  grid  points  respectively  in  the  x,y,z  directions,  which 
were  distributed  over  254  processors  for  parallel  computing.  The  mesh  system  is  more  completely  described 
in  Refs.  [40]  and  [53].  Inflow  profiles  containing  turbulent  structures,  were  generated  from  an  auxiliary  flat- 
plate  simulation  (see  Refs.  [40]  and  [53]  for  details).  The  specification  of  all  boundary  condition  can  also  be 
found  in  Refs.  [40]  and  [53]. 

Active  flow  control,  applied  to  produce  suppression  of  resonant  acoustic  oscillatory  modes  within  the 
cavity,  was  simulated  by  specifying  a  velocity  profile  exiting  through  the  mass-ejection  slot  appearing  Fig.  7. 
This  profile  had  the  following  assumed  functional  form 

v  =  Asin(u;x)  sin2(0.5a^t),  (36) 
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iOx  =TT 


X  —  Xji  \ 

Xj2  —  Xji  J 

cot  =  27 tI  x  5000/uqo. 

In  the  above,  Xji  and  Xj2  are  the  upstream  and  downstream  extents  of  the  slot  respectively,  and  A  is  an 
amplitude  which  could  be  adjusted  to  control  the  mass-flow  rate.  A  dimensional  forcing  frequency  of  5000 
Hz  corresponds  to  the  value  of  tot,  which  matches  that  of  the  experiment.  The  assumed  profile  generates 
a  fluctuating  injection  velocity  which  is  always  positive.  At  the  plane  of  the  jet  exit,  the  pressure  was 
obtained  from  the  inviscid  normal  momentum  equation,  and  the  jet  was  assumed  to  be  isothermal  at  the 
wall  temperature. 

Several  differences  exist  between  the  computed  flowfields,  and  the  experimental  configuration  which  they 
attempt  to  simulate,  observed  in  Fig.  8.  Evident  in  the  figure,  are  sidewall  surfaces  which  mimic  weapons 
bay  doors  in  the  open  position.  The  width  to  length  ratio  of  the  configuration  is  0.2,  which  is  about  twice 
that  of  the  computational  domain.  Mass  injection  was  delivered  by  a  series  powered  resonance  tubes [49] , 
located  beneath  the  jet  exit.  These  tubes  were  fed  by  a  single  plenum  and  discharged  through  the  common 
slot.  The  expelled  flow  was  probably  neither  two-dimensional  nor  isothermal,  but  the  complex  nature  of  the 
interior  region  below  the  slot  was  beyond  the  scope  of  the  simulation. 

5.1  Features  of  the  Time-Mean  Flowfields 

Features  of  the  time-mean  cavity  flowfields  are  described  in  Figs.  9-14.  For  all  of  these  results,  information 
was  obtained  by  spatial  averaging  in  the  homogeneous  direction  (spanwise),  as  well  as  temporally.  Displayed 
in  Fig.  9  are  contours  of  the  mean  pressure  coefficient.  Dark  areas  represent  regions  of  low  pressure,  whereas 
high  pressures  are  lighter.  A  weak  oblique  shock  lies  a  distance  of  approximately  one  cavity  depth  upstream 
of  the  forward  bulkhead  in  the  unsuppressed  case,  and  more  than  twice  that  distance  when  mass  injection 
is  active.  Pressure  levels  within  the  cavity  are  noticeably  lower  (darker)  in  the  suppressed  case.  Found  in 
the  unsuppressed  case  is  a  low  pressure  region  at  the  mouth  of  the  cavity  upstream  of  the  aft  bulkhead. 
For  the  suppressed  case,  the  pressure  in  this  area  is  lower,  and  its  vertical  extent  is  much  greater.  A  more 
quantitative  comparison  of  the  mean  pressure  levels  is  provided  by  the  distributions  along  the  cavity  floor 
(y  =  —0.2)  in  Fig.  10.  Reduction  in  the  level  of  Cp  due  to  suppression  is  apparent. 

Time-mean  streamwise  velocity  contours  appear  in  Fig.  11.  A  thickening  of  the  boundary  layer  upstream 
of  mass  injection  can  be  observed  in  the  suppressed  case.  Also  noted  in  the  suppressed  case,  is  a  reduction 
in  size  of  the  low  speed  region  within  the  cavity  near  the  rear  bulkhead.  Profiles  of  the  mean  streamwise 
velocity  at  three  locations  are  shown  in  Fig.  12.  A  reduction  in  reversed  flow  for  the  suppressed  case  at 
x  =  0.5  is  evident. 

Contours  of  mean  turbulent  kinetic  energy  are  seen  in  Fig.  13.  Although  the  level  of  the  kinetic  energy 
is  elevated  in  the  immediate  vicinity  of  the  injection  slot  due  to  high  frequency  forcing,  its  overall  effect  is  to 
reduce  the  level  downstream  and  within  the  cavity.  Spatial  statistical  characteristics  of  the  cavity  acoustic 
resonance  are  illustrated  by  the  turbulent  kinetic  energy  spanwise  wave-number  spectra  of  Fig.  14.  These 
spectra  were  generated  at  three  streamwise  locations  for  y  =  0.2So/l.  This  position  is  located  within  the 
turbulent  shear  layer  over  the  mouth  of  the  cavity.  The  figure  was  produced  by  constructing  spectra  at 
each  time  step  and  temporally  averaging  the  result.  Here,  the  figure  shows  that  the  energy  is  higher  in  the 
suppressed  case  at  all  streamwise  locations,  most  noticeably  at  x  =  0.2.  It  should  be  noted  that  the  locations 
represented  in  the  figure  lie  within  the  shear  layer  where  forced  injection  was  initiated.  And  although  the 
energy  of  the  shear  layer  increased,  peak  amplitudes  of  the  primary  modes  were  diminished. 

5.2  Features  of  the  Unsteady  Flowfields 

Features  of  the  unsteady  cavity  flowfields  are  elucidated  by  a  series  of  instantaneous  contours  of  flow  variables 
in  Figs.  15  and  17-21.  In  all  of  these  figures,  the  contours  were  formed  at  the  cavity  midspan  centerline. 
Contours  of  the  fluctuating  pressure  p'p'  are  depicted  in  Fig.  15.  A  decrease  in  magnitude  with  suppres¬ 
sion  is  most  noticeable  at  the  rear  bulkhead.  Fluctuating  pressure  profiles  indicated  in  Fig.  16  quantify 
the  mitigating  effect  of  suppression.  Instantaneous  contours  of  Mach  number  in  Fig.  17,  reveal  the  shock 
wave  structures  that  arise  near  the  cavity  rear  bulkhead,  and  upstream  of  the  pulsing  jet.  These  contours 


(37) 

(38) 
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demonstrate  the  shear  layer  that  spans  the  mouth  of  the  cavity.  Dark  contours  in  the  figure  denote  regions 
of  slowly  moving  fluid. 

Presented  in  Figs.  18-21  are  instantaneous  contours  of  the  spanwise  component  of  vorticity.  These 
planar  representations  were  generated  at  the  midspan  location,  for  four  discrete  instants  in  time  ^1,^25^33^4- 
Dynamics  of  the  cavity  fiowfield  are  dominated  by  large-scale  vortical  structures  which  form  aft  of  the 
forward  bulkhead,  and  convect  downstream.  These  structures  evolve  through  a  “roll  up”  of  the  unstable 
shear  layer,  which  is  created  as  the  boundary  layer  leaves  the  surface  ahead  of  the  cavity  at  the  forward  lip. 
The  time  sequence  t\ .  t-2 ,  £3 ,  £4  represents  one  cycle  of  the  vortex  shedding  period,  divided  into  equally  spaced 
increments.  Unsuppressed  and  suppressed  results  were  synchronized,  so  that  at  each  time  instant,  contours 
from  both  cases  corresponded  to  the  identical  point  within  the  vortex  shedding  cycle. 

At  t  =  t\  (Fig.  18),  one  large  vortical  structure  is  forming  in  the  upstream  half  of  the  cavity,  while  another 
has  been  destroyed  as  it  impacted  the  rear  cavity  lip.  Within  the  vortex,  fine-scale  turbulence  is  apparent 
in  the  vorticity  contours.  Fine-scale  turbulence  is  also  observed  in  the  incoming  boundary  layer  upstream 
of  the  forward  lip.  The  figure  shows  that  the  vortex  of  the  suppressed  case  is  weaker  than  its  unsuppressed 
counterpart.  This  is  because  energy  has  been  added  to  the  shear  layer  through  forced  mass  injection,  and 
it  is  better  able  to  withstand  its  natural  tendency  to  roll  up.  In  addition,  the  injection  has  disrupted  the 
coherence  of  the  shear  layer. 

The  state  of  the  fiowfield  at  t  =  t%  is  provided  in  Fig.  19.  Here  the  vortex  has  moved  further  downstream, 
and  now  lies  in  the  center  of  the  cavity.  This  has  produced  a  large  deflection  of  the  shear  layer  in  the 
unsuppressed  case.  The  vortex  begins  its  impact  upon  the  rear  bulkhead  at  t  =  £3  as  is  evident  in  Fig.  20. 
Now,  the  shear  layer  at  the  upstream  cavity  lip  has  returned  to  an  almost  undeflected  condition.  Finally,  at 
t  =  ti  in  Fig.  21,  the  vortex  fully  impinges  upon  the  rear  bulkhead.  The  impingement  is  much  less  severe  in 
the  suppressed  case,  because  the  vortex  is  weaker  and  flatter.  A  new  vortex  can  be  seen  forming  downstream 
of  the  forward  cavity  lip,  thus  completing  the  shedding  cycle.  More  complete  discussions  of  the  instantaneous 
external  shock  system  and  internal  acoustic  wave  propagation  are  given  in  Refs.  [40]  and  [53] . 

Turbulent  kinetic  energy  frequency  spectra  at  three  streamwise  locations  for  y  =  0.2<5q/£  appear  in  Fig.  22. 
The  figure  displays  two  dominant  modes,  particularly  in  the  unsuppressed  case,  where  the  frequency  has 
been  normalized  by  the  first  of  these,  u>\.  A  reduction  in  amplitude  of  the  dominant  modes  is  found  for  all 
streamwise  stations  in  the  suppressed  case.  Addition  of  energy  due  to  mass  injection  appears  for  x  =  0.2, 0.5 
at  oj/w  1  =  27.8  Further  downstream  at  x  =  0.8,  the  effect  of  this  addition  has  diminished.  Similar  to  the 
spatial  wavenumber  spectra  of  Fig.  14,  the  frequency  spectra  were  generated  at  each  spanwise  location,  and 
then  spatially  averaged. 

The  three-dimensional  cavity  fiowfield  is  depicted  in  Fig.  23,  which  provides  instantaneous  total  pressure 
coefficient  contours  at  several  spanwise  stations,  as  well  as  an  iso-surface  of  total  pressure  that  illustrates  the 
primary  vortical  structure.  In  this  representation,  the  spanwise  (z)  direction  has  been  stretched  in  order  to 
more  easily  view  details  of  the  fiowfield.  Although  many  aspects  of  the  flow  have  a  dominant  two-dimensional 
appearance,  the  fine-scale  features  are  clearly  three-dimensional. 

5.3  Comparison  with  Experimental  Data 

Comparison  with  the  experimental  data  of  Refs.  [49,  50,  51,  52]  is  made  in  Figs.  24-27.  This  comparison 
consists  of  fluctuating  pressure  frequency  spectra  at  discrete  points  located  on  the  rear  bulkhead  and  the 
cavity  floor.  As  was  done  for  the  turbulent  kinetic  energy,  frequency  spectra  were  constructed  at  every  z 
location  and  the  results  were  then  spanwise  averaged.  The  amplitude  of  the  fluctuating  pressure  is  presented 
as  sound  pressure  level  ( SPL )  in  dB,  and  the  frequency  ui  is  given  in  Hz,  as  is  customary  for  acoustic 
investigations,  and  is  compatible  with  the  experimental  data.  In  terms  of  the  nondimensional  fluctuating 
pressure,  the  SPL  was  obtained  as 


SPL  =  201og10 


(39) 


where  q  =  2  x  10-5  Pa. 

The  sound  pressure  level  on  the  rear  bulkhead  (x  =  1.0)  at  y  =  —0.04  is  seen  in  Fig.  24.  This  location  is 
just  below  the  aft  lip  of  the  cavity,  and  experiences  the  highest  pressures  due  to  the  oblique  shock  wave  which 
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moves  downstream  preceding  the  vortex,  and  passes  over  the  lip.  For  the  unsuppressed  case,  amplitudes 
of  the  computed  SPL  of  the  two  dominant  modes  compare  well  to  the  experiment.  The  corresponding 
frequencies  of  these  modes  from  the  LES  are  somewhat  smaller  than  those  of  the  measurements.  This  is 
probably  because  the  Reynolds  number  of  the  computations  is  considerably  lower.  It  should  be  noted  that  the 
experimental  frequencies  agree  with  the  empirical  relationship  of  Rossiter[54,  55],  which  is  a  high  Reynolds 
number  correlation.  When  acoustic  suppression  is  active,  both  the  experiment  and  computation  indicate  a 
15  dB  reduction  in  the  amplitude  of  the  dominant  mode. 

Figures  25-27  illustrate  pressure  levels  at  several  streamwise  locations  along  the  floor  of  the  cavity  (y  = 
—0.2).  Amplitudes  on  the  floor  are  not  as  high  as  those  on  the  rear  bulkhead.  Reduction  of  the  acoustic 
modes  with  suppression  is  apparent,  and  the  large-eddy  simulations  compares  favorably  with  the  experiment. 

The  mass  flux  of  forced  injection  for  the  suppressed  case  of  the  experimental  configuration  is  approxi¬ 
mately  145.6  kg/s-m2.  Smaller  values  of  mass  flux  were  also  considered  in  the  investigation,  but  were  not 
effective  in  suppressing  the  acoustic  resonant  modes.  In  the  large-eddy  simulation,  the  time-mean  mass 
flux  was  only  22%  of  the  experimental  value.  When  the  full  experimental  value  of  mass  flux  was  employed 
in  the  computation,  massive  separation  occurred  upstream  of  the  injection  slot.  This  resulted  in  a  large 
displacement  of  the  boundary  layer,  and  the  formation  of  a  near  normal  shock  wave  ahead  of  the  displaced 
region.  The  normal  shock  wave  then  caused  further  separation  of  the  boundary  layer,  and  it  began  to  travel 
upstream.  Eventually  the  shock  reached  the  upstream  boundary,  and  the  computation  was  terminated.  The 
reason  for  this  physical  behavior  was  because  the  Reynolds  number  (Reso)  of  the  computed  flow  was  an  order 
of  magnitude  lower  than  that  of  the  experiment  (see  Table  1).  Therefore,  the  boundary  layer  contained  less 
energy  and  was  not  able  to  withstand  the  disruption  of  strong  mass  injection. 

The  second  mode  of  the  frequency  spectra  found  in  Figs.  25-27  which  occurs  at  approximately  420  Hz 
for  the  LES,  correlates  with  the  shedding  frequency  of  the  vortical  structure.  Because  an  oblique  shock  wave 
is  formed  ahead  of  each  vortex  as  it  is  shed,  the  second  mode  also  corresponds  the  frequency  with  which 
the  shock  impinges  upon  the  rear  bulkhead.  In  addition  however,  there  is  a  system  of  unsteady  pressure 
waves  which  characterize  the  acoustic  resonance  of  the  cavity.  As  each  vortex  forms,  there  eventually  arises 
a  pressure  wave  immediately  beneath  the  vortex  core,  which  travels  downstream  at  about  the  convective 
speed  of  the  vortex.  When  this  wave  impacts  the  aft  bulkhead,  it  is  reflected  from  the  wall  and  then  travels 
upstream.  As  the  wave  travels  upstream,  it  eventually  passes  another  downstream  moving  wave  beneath  the 
next  vortex  which  has  formed.  Eventually,  the  upstream  moving  wave  impacts  the  forward  bulkhead  and  is 
reflected  as  a  downstream  moving  wave.  By  this  time,  yet  another  new  vortex  has  already  begun  to  form 
downstream  of  the  forward  lip.  The  newly  reflect  wave  then  catches  up  with  the  vortex,  and  is  synchronized 
with  its  downstream  travel.  Because  this  entire  cycle  corresponds  to  the  shedding  of  two  vortices,  its 
frequency  is  one  half  that  of  the  second  mode  and  thus  represents  the  first  mode  seen  in  the  spectra.  This 
description,  based  upon  our  observations,  is  consistent  with  those  originally  given  by  Rossiter[54,  55],  but 
differs  in  some  specific  details. 
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6  Leading-Edge  Vortex  Control  on  a  Delta  Wing 

Control  of  the  leading-edge  vortex  generated  about  a  delta  wing  is  a  subject  of  both  fundamental  and 
practical  interest.  Many  swept  wing  flow  control  studies,  utilizing  a  variety  of  actuators  and  strategies,  have 
been  documented  in  the  the  technical  literature  (see  Refs.  [56]  and  [57]).  The  use  of  DBD  plasma  actuators 
offers  an  alternative  approach  for  the  control  of  delta  wing  aerodynamics,  which  requires  no  moving  parts 
or  complex  air  delivery  systems.  At  moderate  Reynolds  numbers,  recent  progress  in  the  design  of  efficient 
plasma  actuators  has  shown  the  potential  for  effective  control  of  such  vortex  flows.  In  the  work  reported 
here,  we  describe  an  investigation  that  explored  control  of  the  complex  flowfield  above  a  slender  delta  wing 
at  high  angle  of  attack  (AOA),  employing  the  empirical  plasma  model  to  represent  plasma-induced  forces. 

The  flow  about  a  flat-plate  delta  wing  with  the  leading  edge  swept  at  75  deg  was  simulated  numer¬ 
ically  at  a  freestream  Mach  number  of  0.1.  Computations  were  performed  for  the  AOA=34  and  38  deg 
and  i?e=25,000,  where  the  characteristic  length  l  was  taken  as  the  root  chord.  For  the  AOA=38  deg,  a 
calculation  was  also  carried  out  at  i?e=9,200  in  order  examine  the  effect  of  Reynolds  number  on  the  location 
of  vortex  breakdown.  The  grid  system  having  an  H-H  topology,  consisted  of  247  x  143  x  209  points  in 
the  streamwise,  spanwise,  and  normal  directions  respectively.  Solutions  were  obtained  with  a  sixth-order 
compact  finite-difference  discretization,  used  in  conjunction  with  an  eighth-order  spatial  filter.  In  order  to 
reduce  computational  resources,  the  flows  were  assumed  to  be  symmetric  about  the  wing  midspan,  so  that 
a  one  half  span  configuration  was  simulated,  with  a  symmetry  plane  along  the  centerline.  A  computational 
time  step  of  At  =  6.25  x  10-5  was  specified  in  order  to  resolve  high  frequency  shear-layer  instabilities  near 
the  wing  apex[58]. 

A  steady  (continuous)  plasma  actuator  with  the  body  force  directed  in  the  downstream  direction,  was 
located  on  the  wing  upper  surface  near  the  apex  as  depicted  schematically  in  Fig.  28.  This  location  was 
found  to  be  most  effective  among  those  which  were  explored,  including  placement  along  the  wing  leading 
and  in  the  trailing  edge  region.  As  noted  previously,  nondimensional  actuator  parameters  were  selected  as 
BC  =  0.018  and  AC  =  0.024  (see  Fig.  6).  The  actuator  was  placed  at  x  =  0.08,  and  spanned  the  region 
0.004  <  y  <  0.024.  The  plasma  scale  parameter  Dc  was  set  to  2400,  which  was  the  same  value  used  previously 
by  Gaitonde  et  al. [33]  for  plasma  control  of  wing  sections. 

At  high  angles  of  attack,  the  baseline  (no  control)  delta  wing  flowfield  is  characterized  by  the  presence  of 
several  interesting  vortical  features[56,  57,  58,  59].  These  include:  1)  boundary-layer  separation  at  the  sharp 
leading  edge,  2)  the  formation  of  a  swept  shear  layer  and  its  roll-up  into  the  primary  leading-edge  vortex, 
and  3)  vortex  breakdown  induced  by  the  adverse  pressure  gradient  near  the  trailing  edge.  Interaction  of 
the  primary  vortex  with  the  boundary  layer  on  the  wing  upper  surface,  gives  rise  to  secondary  vortices.  In 
addition,  the  three-dimensional  “sheet”  feeding  the  leading-edge  vortex  is  known  to  support  both  steady  and 
unsteady  coherent  substructures [58,  60,  61,  62]  associated  with  shear-layer  instabilities,  and  with  unsteady 
boundary-layer  separation  and  vorticity  ejection  on  the  wing  upper  surface. 

6.1  Results  for  AOA=34  deg 

The  body  force  of  the  simulated  DBD  actuator  has  a  profound  impact  on  the  vortex  flow  structure  above 
the  delta  wing.  As  indicated  by  the  time  history  in  Fig.  29,  the  breakdown  location  moves  downstream 
following  the  onset  of  actuation.  After  a  period  of  approximately  three  nondimensional  time  units,  it  has 
traveled  beyond  the  trailing  edge  of  the  wing. 

Instantaneous  flowfields  for  the  baseline  and  control  cases  are  presented  in  Fig.  30  as  contours  in  a 
longitudinal  plane  passing  through  the  center  of  the  vortex.  These  include  pressure  coefficient,  streamwise 
velocity  (u),  and  the  streamwise  component  of  vorticity.  In  the  baseline  case,  the  vortex  core  is  observed  to 
burst  above  the  wing  surface  (Fig.  30a-c).  This  process  is  characterized  by  a  loss  of  both  coherence  and  low 
pressure  in  the  vortex  core,  and  by  the  presence  of  a  significant  region  of  reversed  streamwise  flow.  Although 
the  DBD  actuator  is  placed  within  the  secondary  vortex  region,  it  has  a  significant  effect  on  the  primary 
leading-edge  vortex.  With  actuation  on,  the  vortex  core  remains  intact  over  the  the  entire  length  of  the  wing 
surface  (Fig.  30d-f).  In  addition,  the  core  exhibits  higher  streamwise  velocities  relative  to  the  baseline  case, 
even  prior  to  the  breakdown  location. 

A  perspective  view  of  the  two  flowfields  is  represented  by  instantaneous  cross  plane  contours  of  stream- 
wise  vorticity  and  streamwise  velocity  in  Fig.  31,  and  reveals  well  defined  substructures  and  high  values 
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of  streamwise  vorticity  in  the  shear  layer  of  the  control  case  (Fig.  31a).  It  should  be  noted  in  this  figure 
that  the  left  and  right  halves  of  the  wing  solutions  correspond  to  the  baseline  and  control  simulations  re¬ 
spectively.  The  shear  layer  substructures  are  also  apparent  along  and  above  the  vortex  core  in  Fig.  30f. 
Structural  changes  in  the  shear  layer  induced  by  the  actuator  are  shown  more  clearly  in  Fig.  32,  which 
provides  iso-surfaces  of  instantaneous  streamwise  vorticity.  In  the  baseline  case  (Fig.  32a, b),  unsteady  sub¬ 
structures  are  observed  being  shed  from  the  leading  edge  just  upstream  of  the  breakdown  location.  They 
convect  around  the  leading-edge  vortex,  and  eventually  lose  their  coherence  through  the  interaction  with 
the  expanding  vortex  breakdown  region.  Under  the  influence  of  the  body  force  due  to  actuation,  the  shear 
layer  rearranges  itself  into  a  series  of  steady  helical  substructures  which  emerge  downstream  of  the  actuator 
location  (see  Fig  32d,e).  As  previously  noted,  these  substructures  are  characterized  by  high  values  of  vor¬ 
ticity  that  are  comparable  to  those  in  the  vortex  core.  The  origin  of  steady  substructures  is  linked  to  the 
increase  in  streamwise  velocity  within  the  secondary  vortex,  generated  by  the  actuator.  Figure  32b, e  show 
a  clear  difference  in  the  initial  orientation  of  the  shear  layer  substructures.  The  stationary  structures  are 
inclined  towards  the  wing  trailing  edge  in  the  direction  of  the  prevailing  helical  flow,  whereas  the  axes  of  the 
unsteady  substructures  are  inclined  towards  the  wing  apex. 

A  comparison  of  the  time-mean  flowfields  for  the  baseline  and  control  cases  appears  in  Fig.  33.  On  the 
vertical  plane  through  the  vortex  center  (Fig.  33a),  the  baseline  flow  exhibits  reversed  streamwise  flow,  which 
is  approximately  conical,  over  a  significant  portion  of  the  region  above  the  wing  surface.  With  actuation 
(Fig.  33d),  there  is  a  reduction  of  the  core  streamwise  velocity  just  upstream  of  the  wing  trailing  edge,  but 
the  flow  retains  jet-like  characteristics.  Contours  of  the  root-mean-square  streamwise  velocity  fluctuations 
are  displayed  in  Fig.  33b, e.  High  levels  of  fluctuations  are  present  within  the  vortex  breakdown  region  for 
the  baseline  flow.  Although,  the  vortex  does  not  experience  breakdown  in  the  control  case,  there  is  a  region 
of  significant  fluctuations  near  the  trailing  edge  that  is  associated  with  vortex  core  unsteadiness.  Time- 
mean  shear-layer  structures  are  indicated  in  Fig.  33c, f.  For  the  baseline  case,  these  structures  have  a  weak 
helical  secondary  pattern,  resulting  from  the  secondary  breakdown  of  the  unsteady  vortices,  after  they  are 
shed  from  the  wing  leading  edge.  These  weak  helical  substructures  correspond  in  the  mean,  to  the  imprint 
left  by  vorticity  concentrations  that  were  convected  by  the  prevailing  spiral  flow.  Since  the  control  case 
is  characterized  by  strong  stationary  structures,  the  time-mean  (Fig.  33f)  and  instantaneous  iso-surfaces 
(Fig.  32f)  are  very  similar,  except  near  the  trailing  where  flow  unsteadiness  is  present. 

6.2  Results  for  AOA=38  deg 

At  a  higher  AOA=38  deg,  vortex  breakdown  in  the  baseline  solution  is  closer  to  the  apex  of  the  wing 
(xb  ~  0.36)  as  illustrated  in  Fig.  34a, b.  Application  of  the  steady  plasma  force  displaces  the  breakdown 
downstream  to  Xb  ~  0.6  (Fig.  34c, d).  A  comparison  of  iso-surfaces  of  streamwise  vorticity  in  Fig.  35a, c 
demonstrates  the  transformation  of  the  shear  layer  substructures  under  influence  of  the  DBD  actuator. 
Simulations  at  AOA=38  deg  were  also  carried  out  for  I?e=9,200.  Instantaneous  iso-surfaces  of  streamwise 
vorticity  of  that  computation  are  presented  in  Fig.  35b, d.  An  increase  of  the  relative  magnitude  of  viscous 
forces  results  in  a  smaller  change  of  the  vortex  breakdown  position.  One  intriguing  aspect  of  this  simulation 
is  the  alteration  of  the  breakdown  structure  from  a  spiral  to  a  bubble  type.  For  delta  wing  flows,  this 
situation  is  typically  accompanied  by  a  displacement  of  the  vortex  breakdown  location  towards  the  apex, 
which  is  opposite  of  the  present  situation.  The  abrupt  nature  of  vortex  breakdown  in  the  control  case  is 
also  evident  for  the  higher  Reynolds  number  simulation  in  Fig.  33c.  The  ability  of  plasma-based  actuation 
to  significantly  modify  the  leading-edge  vortex  structure  and  breakdown  location  may  open  new  avenues 
for  enhanced  flight  control.  For  instance,  roll  control  may  be  achieved  without  mechanical  apparatus  by 
activating  DBD  devices  in  an  asymmetric  fashion. 
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7  Efficiency  Enhancement  for  Highly  Loaded  Low-Pressure  Tur¬ 
bine  Blades 

Low-pressure  turbines  are  commonly  utilized  in  the  propulsion  systems  of  unmanned  air  vehicles  employed 
for  reconnaissance  and  combat  purposes.  Due  to  a  decrease  in  atmospheric  density  and  reduced  engine 
speeds  during  high-altitude  cruise,  such  low-pressure  turbines  may  encounter  Reynolds  numbers,  based  upon 
blade  axial  chord  and  inlet  conditions,  below  25,000.  In  this  situation,  boundary  layers  along  a  large  extent 
of  blade  surfaces  can  remain  laminar,  even  in  the  presence  of  elevated  freestream  turbulence  levels.  These 
laminar  boundary  layers  are  then  particularly  susceptible  to  flow  separation  over  the  aft  portion  of  blade 
suction  surfaces,  promoted  by  the  adverse  pressure  gradient  as  the  flow  turns.  This  results  in  a  breakdown  of 
the  flow,  transition  to  turbulence,  blockage  in  flow  passages,  and  a  significant  reduction  in  turbine  efficiency. 

A  number  of  experimental  investigations  by  Bons  et  al.  [63,  64,  65]  and  by  Sondergaarcl  et  al.[66,  67] 
have  explored  the  use  of  both  steady  and  pulsed  vortex  generator  jets,  which  may  be  actuated  upon  demand, 
as  a  means  of  flow  control  in  low-pressure  turbines.  In  particular,  the  work  by  Sondergaarcl  et  al. [67] 
considered  the  feasibility  of  increasing  the  inter-blade  spacing,  thereby  raising  the  per  blade  loading.  For 
practical  applications,  a  higher  loading  can  reduce  the  turbine  part  count  and  stage  weight.  Increased 
blade  spacing  however,  is  accompanied  by  more  extensive  boundary-layer  separation  on  each  blade  due  to 
uncovered  turning,  resulting  in  a  further  reduction  of  efficiency  and  additional  wake  losses.  More  highly 
loaded  configurations  are  therefore  only  functional  when  used  in  conjunction  with  flow  control.  Rizzetta 
and  Visbal[68,  69,  36,  70,  71]  carried  out  a  series  of  numerical  simulations  employing  both  vortex  generator 
jets  and  plasma  actuators  to  mitigate  separation  and  improve  efficiency  for  a  highly  loaded  linear  cascade 
of  low-pressure  turbine  blades,  that  is  similar  to  the  experiments  cited  above.  Either  of  these  devices 
may  be  activated  upon  demand,  remaining  idle  at  sea  level  where  flow  along  the  blades  is  attached,  and 
being  actuated  at  altitude  when  separation  occurs.  Complete  details  of  the  computations  are  found  in 
Refs.  [68,  69,  36,  70,  71],  and  some  of  the  significant  results  are  summarized  below. 

Shown  in  Fig.  36  is  a  schematic  representation  of  the  turbine  blade  shape,  given  by  the  Pratt  &  Whitney 
“PakB”  research  design,  which  is  a  Mach  number  scaled  version  of  geometries  typically  used  in  highly  loaded 
low-pressure  turbines. [63,  64,  65,  67,  66]  This  blade  geometry  has  an  inlet  flow  angle  a*  =  35.0  deg,  and 
a  design  exit  flow  angle  a0  =  60.0  deg.  The  axial  chord  to  spacing  ratio  (solidity)  is  0.75,  resulting  in  an 
inter-blade  spacing  B  =  4/3. 

To  conserve  computational  resources,  only  a  single  turbine  blade  passage  was  considered,  and  periodic 
conditions  were  enforced  in  the  vertical  direction  (y)  to  represent  a  single  turbine  stage.  A  mesh  system 
consisting  of  several  distinct  overlapping  grids  was  used  to  define  the  flowfield,  which  are  observed  in  Fig.  37. 
Figure  37a  depicts  the  O-gricl  topology  of  the  basic  grid  about  the  turbine  blade,  which  was  comprised  of  348 
points  in  the  circumferential  direction  (/),  189  points  in  the  blade-normal  direction  ( J),  and  101  points  in  the 
spanwise  direction  (K).  Periodic  boundaries  along  I\u  —  l2u  and  In  — 121  are  apparent  in  the  figure.  Periodic 
conditions  were  also  enforced  in  the  spanwise  direction.  An  embedded  locally  refined  region  (Fig.  37b)  of 
(313  x  185  x  101)  grid  points  in  (/,  J ,  K)  respectively,  was  used  to  capture  the  correct  fluid  physics  in  the 
actuator  and  near-wall  regions.  In  order  to  facilitate  application  of  inflow  and  outflow  conditions  to  the 
turbine  blade  domain,  overset  grids  were  utilized  upstream  and  downstream  of  the  blade  region,  and  are 
provided  in  Fig.  37c.  The  total  number  of  grid  points  in  all  domains  was  approximately  12.0  x  106. 

The  planar  grids  appearing  in  Fig.  37  were  uniformly  distributed  across  the  span,  where  the  nondimen- 
sional  spanwise  extent  of  the  computational  domain  is  s.  In  the  case  of  vortex-generating  jets,  s  was  taken 
as  the  inter-jet  spacing  of  0.112.  For  the  baseline  (no  control)  and  plasma-based  control  simulations,  s=0.2. 
This  value  of  s  was  found  to  be  adequate  to  capture  the  transitional  flowfield  in  the  prior  investigation  of 
Ref.  [37].  Flow  variables  in  the  overlapping  mesh  regions  of  Fig.  37  were  obtained  from  the  aforementioned 
explicit  Lagrangian  interpolation  formulae,  with  sixth-order  spatial  accuracy.  The  turbine  blade  axial  chord 
was  selected  as  the  characteristic  reference  length,  and  the  reference  Mach  number  was  0.1.  Blockage  in  the 
flow  passage  caused  the  flow  to  be  retarded,  lowering  the  Mach  number  at  the  inflow  boundary.  The  reference 
Reynolds  was  then  iteratively  determined  in  order  to  match  conditions  at  the  inflow,  which  corresponded  to 
previous  experiments  and  computations.  A  more  thorough  explanation  of  this  procedure,  and  details  of  all 
boundary  conditions,  can  be  found  in  Refs.  [68,  69,  36,  70,  71].  The  Reynolds  number  based  upon  inflow 
conditions  was  approximately  25,000. 

All  solutions  were  obtained  using  a  fourth-order  compact  differencing  scheme  and  a  sixth-order  spatial 
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filter,  with  a  nondimensional  time  step  At  =  1.5  x  10~4.  Unsteady  actuation,  either  vortex-generating  jets 
or  plasma  control,  consisted  of  a  period  of  time  when  devices  were  active,  followed  by  an  interval  when  they 
were  idle.  The  duty  cycle  represents  the  portion  of  the  total  duration  over  which  control  is  active,  and  was 
50%  for  all  results  presented  here.  The  analytical  description  for  an  amplitude  function  employing  the  50% 
duty  cycle  consisted  of  a  series  of  piecewise  continuous  cubic  and  linear  functions,  and  is  displayed  in  Fig.  38. 
Here,  td  is  the  portion  of  the  fundamental  period  tp  over  which  the  device  is  active.  The  duty  cycle  is  given 
by  the  ratio  td/tp  x  100  expressed  as  a  percentage.  This  amplitude  function  is  then  appended  as  a  factor  to 
the  vortex  jet  exit  velocity  or  plasma  force  magnitude  in  order  to  create  pulsed  control.  It  should  be  noted 
that  the  applied  waveform  introduces  multiple  harmonics  of  the  primary  frequency,  as  was  demonstrated  in 
Refs.  [72]  and  [35].  The  period  illustrated  in  Fig.  38  was  represented  by  1300  time  steps,  corresponding  to  a 
nondimensional  frequency  of  7.56. 

From  computed  time-mean  flowfields,  control  effectiveness  was  quantified  by  calculation  of  the  integrated 
wake  total  pressure  loss  coefficient  Cw,  defined  as 
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Equation  40  was  evaluated  along  the  upstream  boundary  of  the  downstream  mesh  in  Fig.  37c,  which  is 
located  0.67  chords  downstream  of  the  blade  trailing  edge. 


7.1  Results  for  Vortex-Generating  Jets 

In  the  above-mentioned  experimental  studies[63,  64,  65,  66,  67],  jets  were  created  by  blowing  air  through 
holes  which  had  been  drilled  in  the  blade  surface  at  a  pitch  angle  of  30  deg  and  a  skew  angle  of  90  deg.  The 
jets  were  situated  at  the  time-mean  separation  location,  and  the  geometry  is  indicated  in  Fig.  39.  Here,  the 
pitch  is  defined  as  the  angle  the  jet  makes  with  the  local  surface,  and  the  skew  is  the  angle  of  the  projection  of 
the  jet  on  the  surface,  relative  to  the  “local  freestream”  direction.  [65]  The  size  of  the  drill  used  to  develop  the 
holes  is  commonly  referred  to  as  the  jet  diameter,  which  was  0.001  m.  Because  of  the  orientation  however, 
the  jet  exit  geometric  shape  is  elliptic  as  seen  in  Fig.  39,  and  the  jet  exit  velocity  vector  has  components 
only  in  the  blade-normal  and  spanwise  (z)  directions.  In  order  to  simulate  the  flow  within  the  jet  nozzle, 
additional  overset  meshes  were  employed  in  the  jet  exit  region,  and  interior  to  the  blade  geometry.  These 
are  evident  in  Fig.  39.  A  description  of  inflow  conditions  to  the  jet  nozzle  is  given  in  Refs.  [68]  and  [69]. 

Time-mean  surface  pressure  coefficient  distributions  for  the  baseline  (no  control)  and  flow  control  compu¬ 
tations  are  presented  in  Fig.  40.  For  the  flow  control  cases,  these  distributions  were  obtained  at  the  spanwise 
location  of  the  periodic  boundary  between  two  jets.  As  indicated  in  the  figure,  the  jet  was  positioned  at 
x  =  0.37,  which  coincides  with  the  separation  point  of  the  baseline  solution.  The  large  plateau  region  in  the 
baseline  distribution  is  characteristic  of  massively  separated  flow.  Because  of  reduced  blockage  and  increased 
inflow  velocity,  the  effect  of  flow  control  is  to  decrease  the  pressure  on  the  upstream  portion  of  the  suction 
surface,  while  increasing  it  downstream,  relative  to  the  baseline  case.  Separation  occurs  at  x  =  0.56  for  the 
steady  blowing  result,  and  at  x  =  0.58  for  pulsed  injection.  Only  minor  differences,  most  noticeable  near  the 
jet,  are  observed  between  the  two  flow  control  solutions.  Also  seen  in  the  figure,  is  a  coarse-mesh  distribution 
for  the  baseline  case,  which  indicates  little  sensitivity  to  grid  resolution,  except  near  the  trailing-edge  region. 
More  extensive  results  of  the  grid  resolution  study  for  the  baseline  case  are  found  in  Ref.  [73]. 

Time-mean  contours  of  the  streamwise  component  of  velocity  and  spanwise  component  of  vorticity  are 
shown  in  Fig.  41.  Contours  for  the  flow  control  results  were  taken  at  the  plane  of  the  periodic  boundary 
between  jets,  while  those  of  the  baseline  case  have  been  spanwise  averaged.  An  increased  region  of  attached 
flow  due  to  flow  control  is  apparent,  resulting  in  a  decrease  of  the  wake  thickness.  Maintaining  attached  flow 
and  decreasing  the  extent  of  the  wake  are  exhibited  by  the  flow  control  results. 

Spanwise  turbulent  kinetic  energy  wave-number  spectra  are  provided  in  Fig.  42.  These  spectra  were 
generated  along  lines  in  the  z  direction  at  a  nondimensional  distance  of  n  =  0.03  from  the  blade  surface. 
This  location  ( n  =  0.03)  is  approximately  equal  to  one  half  of  the  boundary-layer  thickness  of  the  time-mean 
velocity  profile  upstream  of  separation.  At  the  two  most  upstream  stations  (x  =  0.50, 0.70),  Ekz  is  higher  for 
the  flow  control  cases  than  that  of  the  baseline  due  to  energy  being  added  to  the  flow.  In  the  downstream 
region  (x  =  0.90),  flow  control  has  mitigated  separation  and  the  associated  breakdown  into  a  more  chaotic 
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situation,  so  that  the  turbulent  kinetic  energy  of  the  baseline  case  is  higher.  Because  of  the  low  Reynolds 
number,  only  a  small  portion  of  the  spectrum  at  low  wave  number  lies  in  the  inertial  range. 

Instantaneous  contours  of  the  streamwise  velocity  component  and  spanwise  vorticity  appear  in  Fig.  43. 
Once  again,  the  baseline  result  was  extracted  at  the  midspan  location,  while  those  of  the  flow  control  cases 
were  situated  between  control  jets.  Similar  to  the  time-mean  contours  of  Fig.  41,  it  is  noted  that  flow  control 
maintains  attached  flow  and  decreases  the  vertical  extent  of  the  wake  relative  to  the  baseline  case.  Unsteady 
structures  are  visible  in  separated  flow  regions.  It  is  seen  in  the  flow  control  cases  that  vorticity  (Fig.  43b) 
is  being  generated  in  the  boundary  layer  in  close  proximity  to  the  blade  surface.  When  the  boundary  layers 
separate  in  these  situations,  it  is  much  less  dramatic  than  in  the  uncontrolled  case.  For  the  baseline  solution, 
the  extensive  unsteady  separated  flow  region  has  a  richer  content  of  small  scale  structures  due  to  breakdown 
and  a  transition  to  a  more  complex  situation. 

Turbulent  kinetic  energy  frequency  spectra  are  observed  in  Fig.  44.  The  data  used  to  generate  these 
spectra  was  collected  at  n  =  0.03,  similar  to  that  of  Fig.  42.  Like  the  spanwise  wave  number  spectra  of 
Fig.  42,  the  frequency  spectra  for  the  flow  control  cases  is  substantially  higher  than  that  of  the  baseline 
result  at  the  upstream  stations  (x  =  0.50,  0.70)  because  of  energy  addition.  At  x  =  0.50,  the  discrete  peaks 
in  Ew  for  the  pulsed  injection  case  correspond  to  harmonics  of  the  forcing  frequency.  The  occurrence  of 
this  behavior  was  described  in  detail  in  Ref.  [73] .  As  noted  previously,  turbulent  energy  of  the  baseline  case 
exceeds  that  of  the  flow  control  results  near  the  blade  trailing  edge  (x  =  0.90). 

A  three-dimensional  representation  of  the  flow  is  depicted  by  iso-surfaces  of  vorticity  magnitude  in  the 
trailing-edge  region  seen  in  Fig.  45.  The  value  of  iso-surfaces  correspond  approximately  to  that  at  the  edge 
of  the  shear  layer  upstream  of  separation.  Both  the  vertical  and  spanwise  extent  of  the  turbulent  structures 
are  visible  for  each  case. 

For  both  the  steady  and  pulsed  cases,  the  fundamental  effect  of  vortex-generating  jets  was  to  energize 
the  blade  boundary  layer  due  to  the  transfer  of  fluid  momentum  and  mixing.  This  helped  maintain  attached 
flow  along  the  blade  surface  for  a  distance  of  19-21%  greater  than  that  of  the  baseline  case.  As  a  result,  wake 
total  pressure  losses  were  decreased  by  53-56%.  In  the  pulsed  case,  mixing  was  enhanced  through  unsteady 
forcing.  Although  no  forcing  was  employed  with  the  steady  case,  the  jet  was  inherently  unstable,  and  broke 
down  shortly  after  exiting  the  nozzle.  This  behavior  also  served  to  increase  mixing.  Because  of  the  50% 
duty  cycle,  the  pulsed  case  required  less  mass  flow  than  the  steady  jet  to  achieve  a  similar  improvement  in 
efficiency. 

7.2  Results  for  Plasma-Based  Control 

Extensive  simulations  for  the  highly  loaded  low-pressure  turbine  blade  were  performed  by  Rizzetta  and 
Visbal[36,  70,  71]  in  order  to  explore  plasma-based  flow  control  strategies.  During  the  study,  many  aspects 
of  control  were  investigated  including  the  magnitude  of  the  plasma  force  required  to  establish  control,  the 
use  of  both  continuous  and  pulsed  actuation,  the  effect  of  forcing  frequency  and  duty  cycle  for  pulsing,  the 
location  for  the  placement  of  actuators,  the  use  of  both  full  span  and  finite  distributed  arrays  of  actuators, 
and  the  direction  of  the  plasma-induced  force.  For  the  later  of  these,  the  actuators  were  oriented  such  that 
the  force  was  directed  in  the  coflowing,  counterflowing,  or  spanwise  directions.  A  primary  objective  of  these 
studies  was  to  determine  effective  plasma-control  arrangements  which  required  minimal  power  usage.  Some 
results  of  these  efforts  are  reported  below,  all  of  which  were  obtained  for  pulsed  actuation  with  a  value  of 
the  plasma  scale  parameter  Dc  =  25,  that  represents  a  low  energy  demand. 

Configurations  considered  in  these  results  are  represented  in  Fig.  46,  where  the  arrows  indicate  the 
actuation  direction.  It  can  be  noted  in  the  figure,  that  cases  A,  B,  C,  and  D  correspond  to  counterflow 
actuation,  while  cases  E  and  F  have  spanwise  oriented  actuators.  Because  configurations  E  and  F  were  less 
effective  arrangements,  no  results  for  those  cases  will  be  presented  here. 

Time-mean  surface  pressure  coefficient  distributions  are  displayed  in  Fig.  47.  These  results  were  obtained 
along  the  centerline  of  the  computational  domain  for  each  solution.  The  large  plateau  region  in  the  baseline 
distribution  is  characteristic  of  a  massively  separated  flow,  as  was  noted  in  Fig.  40  for  the  vortex-generating 
jet  solutions.  Again,  the  dominant  effect  of  flow  control  is  to  reduce  separation  of  the  time-mean  flowfield 
on  the  blade  suction  surface.  Cases  A,  B,ancl  C  effectively  maintained  attached  flow  to  the  blade  trailing 
edge,  thus  completely  eliminating  the  plateau.  Although  the  finite-span  case  (configuration  D)  reduced  the 
extent  of  separation,  a  smaller  plateau  is  still  present. 
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A  composite  visual  comparison  of  the  flowfields  for  these  cases  is  exhibited  in  Fig.  48.  Contours  of 
streamwise  velocity  are  displayed  in  the  top  row  of  the  figure,  streamlines  in  the  middle  row,  and  pressure 
coefficient  (Cp)  at  the  bottom.  Massive  separation  in  the  baseline  case  is  apparent.  Results  depicted  here 
are  consistent  with  those  of  the  surface  pressure  distributions  (Fig.  47).  Streamlines  illustrate  the  massively 
separated  flow  of  the  baseline  case.  Completely  attached  time-mean  flow  for  cases  A,  B,  and  C  is  apparent 
in  the  figure.  A  reduction  of  the  separated  flow  region  is  evident  for  case  D.  The  thin  wakes  for  cases  A,  B, 
and  C,  relative  to  the  baseline  flow,  are  visible  in  the  streamwise  velocity  contours.  The  effect  of  separated 
flow  regions  occurring  for  the  baseline  case  and  configuration  D,  is  reflected  in  the  Cp  contours. 

Spanwise  turbulent  kinetic  energy  wave-number  spectra  are  shown  in  Fig.  49.  These  spectra  were  gen¬ 
erated  at  n  =  0.03,  similar  to  those  of  Fig.  42.  Because  the  spectra  for  configuration  B  was  identical  to 
that  of  C,  and  case  D  was  similar  to  the  baseline,  those  results  are  not  presented.  At  the  most  upstream 
station  (a:  =  0.5),  case  C  contains  more  energy  than  either  the  baseline  or  configuration  A  solutions.  This 
is  also  true  further  downstream  at  x  =  0.7.  And  although  energy  is  being  added  to  the  flow,  it  is  believed 
this  situation  reflects  an  early  transition  to  turbulence.  Near  the  trailing  edge  ( x  =  0.9),  all  situations  are 
similar,  toward  the  end  of  the  transition  process.  It  should  be  noted  that  these  comparisons  are  performed 
at  a  location  very  close  to  the  blade  surface.  Away  from  the  blade,  the  total  energy  contained  in  the  large 
structures  of  the  massively  separated  region  for  the  baseline  case  is  very  large.  As  noted  previously,  only 
a  small  portion  of  the  spectrum  at  low  wave  number  for  x  =  0.9  lies  in  the  inertial  range  due  to  the  low 
Reynolds  number. 

Provided  in  Fig.  50  are  instantaneous  streamlines  at  the  blade  midspan  for  configuration  A.  Seen  in  the 
figure,  is  the  flowfield  near  the  actuator  location  at  four  frames  during  the  pulsing  cycle.  For  t/tp  =  0.0, 
control  has  been  inactive  for  one  half  of  the  pulsing  cycle,  and  the  flow  is  attached.  At  t/tp  =  0.2,  actuation 
has  taken  place  for  20%  of  the  cycle,  creating  a  small  separation  bubble  at  the  blade  surface.  By  t/tp  =  0.4, 
the  region  of  reversed  flow  has  grown  in  size.  Finally,  when  t/tp  =  0.6  the  active  portion  of  the  cycle  has 
ended,  and  the  separation  bubble  begins  to  convect  downstream. 

A  composite  representation  of  the  instantaneous  flowfields  for  these  cases  appears  in  Fig.  51.  Contours  of 
streamwise  velocity  and  spanwise  vorticity  along  the  centerline  plane  are  observed  in  the  top  and  middle  rows 
of  the  figure  respectively,  while  spanwise  vorticity  on  the  blade  surface  viewed  from  above,  is  at  the  bottom. 
The  middle  row  of  spanwise  vorticity  is  particularly  useful  for  understanding  the  control  mechanisms.  It 
was  shown  in  Ref.  [36]  for  configuration  A,  that  control  was  imposed  by  modifying  the  inherently  unstable 
boundary  layer  near  its  separation  point.  The  plasma  actuator  forced  the  boundary  layer  to  roll  up  into 
small  vortices  just  downstream  of  the  actuator  location.  These  vortices  were  shed  at  a  frequency  one  half 
that  of  the  pulsed  control,  and  convected  downstream  at  a  distance  from  the  blade  surface  which  was  not 
large.  This  greatly  enhanced  mixing,  and  brought  higher  momentum  fluid  into  the  boundary  layer,  thereby 
maintaining  time-mean  attached  flow  and  reducing  wake  losses.  The  fairly  coherent  vortical  structures  are 
visible  in  the  spanwise  vorticity  contours  (middle  row).  Dark  portions  of  the  surface  vorticity  (bottom  row) 
signify  regions  of  attached  flow. 

For  configuration  B,  actuation  was  applied  further  upstream  from  separation  than  in  configuration  A,  and 
expedited  transition.  As  a  result,  no  coherent  vortices  evolved,  and  the  flowfield  was  dominated  by  fine-scale 
fluid  structures.  Configuration  C  was  specifically  designed  to  eliminate  fluid  coherence.  It  was  believed  that 
coherence  could  promote  structural  fatigue,  and  be  detrimental  to  instrumentation  and  acoustic  signature 
for  turbine  applications.  The  approach  with  configuration  C  was  to  apply  a  second  counterflow  actuator 
downstream  of  the  first,  to  help  break  down  coherent  structures.  In  addition,  the  second  actuator  was  of 
finite  span,  so  as  to  reduce  two-dimensional  effects.  It  is  noted  for  configuration  C  in  Fig.  51  that  only 
fine-scale  structures  are  present.  Both  configurations  B  and  C  have  wake  total  pressure  loss  coefficients  that 
are  slightly  less  than  that  of  configuration  A. 

Configuration  D  considered  only  a  single  finite-span  counterflow  actuator,  at  the  same  location  as  that 
of  configuration  A.  Although  some  control  was  exerted  in  this  situation,  the  actuation  was  less  effective. 
The  arrangement  failed  to  manipulate  the  boundary  layer  in  a  useful  manner,  or  to  generate  transition. 
Effectiveness  arose  only  through  enhanced  mixing. 

Turbulent  kinetic  energy  frequency  spectra  at  n  =  0.03  for  configurations  A  and  B  are  displayed  in 
Fig.  52.  As  occurred  for  pulsed  vortex-jet  control,  discrete  peaks  in  Eu  near  the  actuator  location  ( x  =  0.5) 
with  plasma  actuation,  correspond  to  harmonics  of  the  pulsing  frequency.  At  x  =  0.7,  sub-harmonics  of  the 
pulsing  frequency  have  been  excited  in  the  spectra  of  configuration  A.  This  phenomenon  was  caused  by  the 
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small  vortices  generated  in  the  boundary  layer,  downstream  of  the  actuator.  No  sub-harmonics  are  apparent 
for  configuration  B  (or  configuration  C,  not  shown),  as  the  small  coherent  vortical  structures  have  been 
suppressed.  Turbulence  levels  of  configuration  A  are  higher  than  the  baseline  at  x  =  0.5  and  x  =  0.7  due  to 
energy  being  added  to  the  flowfield.  Peak  levels  of  configuration  B  are  higher  than  those  of  configuration  A 
at  these  location,  although  the  broadband  content  is  lower.  Farther  downstream  in  the  trailing-edge  region 
(x  =  0.9),  all  spectra  attain  a  similar  level.  Peaks  associated  with  actuation  in  configuration  A  at  this 
location,  are  not  present  for  configurations  B  (or  configuration  C,  not  shown).  A  small  inertial  range  in  the 
spectra  can  be  observed  in  all  results. 

A  three-dimensional  representation  of  the  flow  in  the  trailing-edge  region  for  configurations  A  and  C  is 
depicted  in  Fig.  53  by  isosurfaces  of  vorticity  magnitude,  which  have  been  colored  by  the  streamwise  velocity. 
The  value  of  the  isosurfaces  corresponds  approximately  to  that  at  the  edge  of  the  shear  layer  upstream  of 
separation.  Massive  separation  is  evident  in  the  baseline  flow,  as  are  two-dimensional  coherent  structures  for 
configuration  A.  The  elimination  of  these  structures  in  configuration  C  is  noticeable.  Configurations  A,  B, 
and  C  were  able  to  achieve  a  reduction  in  wake  total  pressure  losses  of  81-86%.  This  was  better  than  that 
attained  with  vortex  generating  jets.  For  configuration  D,  losses  were  only  lowered  by  39%. 
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Table  2:  Hump  flow  reattachment  locations 


case 

experiment  [75,  76] 

LES 

RANS 

baseline 

1.11  ±0.003 

1.139 

1.23 

suction 

0.94  ±0.005 

0.978 

1.14 

blowing/suction 

0.98 

1.097 

1.29 

8  Separation  Control  for  a  Wall-Mounted  Hump  Model 

The  control  of  separation  for  wall-bounded  flows  is  an  important  and  challenging  area  for  numerical  simu¬ 
lation.  In  order  to  define  some  of  these  challenges  and  assess  current  capabilities,  the  CFD  Validation  on 
Synthetic  Jets  and  Turbulent  Separation  Control  Workshop,  sponsored  by  NASA  Langley  Research  Center, 
was  held  at  Williamsburg,  VA  in  March  of  2004.  [74]  The  main  objective  of  the  workshop  was  to  develop 
a  comprehensive  database  of  experiments  employing  flow  control  for  subsequent  validation  and  comparison 
with  various  computational  simulation  techniques.  One  of  the  cases  at  this  workshop  corresponded  to  flow 
over  a  wall-mounted  hump  model,  which  was  investigated  experimentally[75,  76]  with  and  without  flow 
control.  In  the  control  cases,  both  steady  suction  and  zero  net  mass  flux  blowing/suction  were  considered. 

The  hump  geometry  simulates  the  upper  surface  of  a  20%  thick  Glauert-Goldschmied  airfoil,  with  a  chord 
of  0.4200  meters,  a  maximum  height  of  0.0537  meters,  and  a  span  of  0.5842  meters.  A  slot  was  situated 
at  the  65%  chord  location,  and  connected  to  an  internal  plenum  for  control  purposes.  The  geometry  of  the 
configuration  and  overset  mesh  system  are  illustrated  in  Fig.  54.  The  hump  portion  of  the  computational 
domain  was  described  by  (818  x  151  x  165)  points  in  the  (/,  J,  K)  directions,  while  a  grid  of  (41  x  133  x  165) 
defined  the  plenum.  Both  steady  suction  and  oscillatory  blowing/suction  flow  control  were  applied  via 
boundary  conditions  enforced  at  the  inlet  to  the  plenum.  Complete  details  of  the  simulations  are  documented 
in  Refs.  [77,  78,  79], 

Experimental  measurements  were  conducted  at  Re  =  9.36  x  105,  but  the  simulations  were  performed 
for  Moo  =  0.10  and  Re  =  2.0  x  105,  where  the  model  chord  was  used  as  the  characteristic  length.  A  lower 
Reynolds  number  was  employed  in  the  computations  to  maintain  LES  resolution  with  acceptable  levels  of 
resource  requirements.  The  large-eddy  simulations  utilized  a  fourth-order  compact  difference  scheme  and  a 
sixth-order  spatial  filter.  Although  not  central  to  the  main  purpose  of  the  computations,  RANS  simulations 
were  also  carried  out  and  are  shown  below  for  comparison.  A  full  description  of  the  RANS  approach  can  be 
found  in  Refs.  [80]  and  Morgan:AIAAJ2006. 

8.1  Features  of  the  Time-Mean  Flowfields 

A  comparison  of  the  time-mean  surface  pressure  coefficients  is  presented  in  Fig.  55.  The  experimental  data, 
LES,  and  RANS  solutions  all  display  a  large  peak  near  the  mid-chord  of  the  hump,  which  is  associated 
with  flow  acceleration.  A  small  plateau  region  downstream  of  the  flow  control  slot  is  characteristic  of  flow 
separation.  All  simulations  compare  well  with  the  experiment  upstream  of  the  control  slot,  where  the  flow  is 
attached.  The  LES  solutions  compare  reasonable  well  with  the  experiment  in  the  separated  flow  and  wake 
regions.  Deficiencies  of  the  RANS  computations  are  particularly  evident  for  the  suction  and  oscillatory  cases. 

Comparisons  of  time-mean  streamwise  velocity  contours  are  indicated  in  Figs.  56-58.  For  the  baseline 
case  (Fig.  56),  the  LES  and  experiment  have  stronger  reverse  flow  than  the  RANS,  in  the  near- wall  region  for 
0.8  <  x  <  1.0.  The  main  difference  between  the  LES  and  the  experiment,  is  that  the  region  of  fastest  reverse 
flow  occurs  5%  of  the  chord  further  upstream  for  the  LES  than  is  observed  in  the  experiment.  The  RANS 
solution  clearly  shows  a  longer  separation  bubble  than  either  the  experimental  or  LES  results.  Reduction 
of  the  separated  flow  region  by  applying  suction  is  apparent  in  Fig.  57.  Control  appears  less  effective  with 
blowing/suction  in  Fig.  58.  For  both  control  cases,  the  RANS  solution  predicts  a  larger  separated  flow 
region  than  either  the  LES  or  the  experiment,  which  compare  favorably.  Locations  of  the  time-mean  flow 
reattachment  points  are  provided  in  Table  2. 
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8.2  Features  of  the  Unsteady  Flowfields 

Contour  plots  comparing  the  time-mean  Reynolds  stresses  are  shown  in  Figs.  59-61.  For  the  baseline  case 
(Fig.  59),  it  can  be  clearly  seen  that  the  LES  developed  a  smaller,  thinner  region  of  minimum  Reynolds 
stress  which  was  located  further  upstream  than  in  the  experiment.  The  RANS  solution  does  not  reach  the 
same  levels  of  Reynolds  stress  as  that  of  the  experiment  and  LES.  The  region  of  minimum  magnitude  of  the 
Reynolds  stresses  for  the  RANS  result  is  significantly  smaller  and  is  located  further  downstream.  Reynolds 
stress  profiles  for  steady  suction  are  compared  in  Fig.  60.  Overall  the  LES  and  experimental  contours  have 
the  same  qualitative  shape  and  magnitude.  The  LES  however,  did  not  develop  the  lowest  magnitude  of 
Reynolds  stress  that  occurred  in  the  experiment.  Contours  of  the  RANS  solution  display  generally  lower 
values  of  Reynolds  stress.  Exhibited  in  Fig.  61  are  contour  plots  of  Reynolds  stress  for  the  blowing/suction 
case.  Here,  the  LES  fails  to  attain  the  same  minimum  levels  of  the  experiment.  This  is  attributed  to 
the  circumstance  that  the  same  magnitudes  of  blowing/suction  were  employed  for  the  LES  as  that  of  the 
experiment,  but  the  Reynolds  number  was  considerably  lower. 

In  order  to  represent  characteristics  of  the  unsteady  flowfield  for  the  oscillatory  blowing/suction  case, 
phase-averaged  planar  contours  of  the  spanwise  component  of  vorticity  were  collected  experimentally.  These 
are  compared  with  the  computations  for  four  values  of  the  phase  angle  in  Figs.  62-65.  The  phase  angle 
0=90  deg  corresponds  to  the  maximum  blowing  rate  and  0=270  deg  to  maximum  suction  in  the  oscillatory 
cycle.  Overall,  both  the  LES  and  RANS  solutions  have  good  agreement  with  the  experiment  in  the  general 
behavior  of  the  flow,  including  size  and  locations  of  the  vortex  generated  by  the  pulsing.  The  RANS  results 
display  more  coherent  structures  throughout  the  blowing  and  suction  cycle,  which  do  not  dissipate  as  rapidly 
in  the  wake. 

Instantaneous  streamwise  velocity  contours  in  the  x  —  y  plane  and  iso-surfaces  of  vorticity  magnitude  for 
the  blowing/suction  case  are  seen  in  Fig.  66.  We  note  that  this  figure  has  been  stretched  by  a  factor  of  two 
in  the  spanwise  direction  to  enhance  visualization.  The  turbulent  boundary  layer  with  multiple  spanwise 
structures  can  be  seen  in  the  upstream  region.  These  turbulent  structures  diminish  as  the  flow  accelerates 
over  the  hump.  After  separation,  the  effect  of  pulsing  is  evident  in  the  large  spanwise  structures  that  convect 
downstream. 
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9  Summary  and  Conclusion 

The  description  of  a  high-order  finite-difference  scheme,  suitable  for  large-eddy  simulations  which  consider 
active  flow  control,  has  been  presented.  The  numerical  method  is  predicated  upon  an  implicit  approximately- 
factored  time-marching  approach,  employing  Newton-like  subiterations,  that  is  efficient  for  wall-bounded 
flows.  Spatial  derivatives  are  evaluated  via  an  implicit  compact  approximation,  thereby  reducing  the  asso¬ 
ciated  dispersion  error.  The  differencing  technique  is  use  in  conjunction  with  a  similar  implicit  high-order 
low-pass  Pade-type  spatial  filter,  that  augments  robustness  by  adding  dissipation  only  for  under-resolved 
high-frequency  spatial  modes  of  the  solution.  Because  of  this  selective  nature  of  the  filtering  process,  it  also 
serves  as  an  implicit  subgrid-scale  stress  model.  Flexibility  of  the  method  for  application  to  geometrically 
complex  configurations  is  provided  by  an  overset  grid  approach.  This  allows  structured  meshes  to  be  utilized 
and  high-order  spatial  accuracy  to  be  achieved.  A  high-order  Lagrangian  interpolation  procedure  main¬ 
tains  accuracy  at  overlapping  grid  boundaries.  Domain  decomposition  is  applied  for  processing  on  parallel 
computing  platforms. 

Several  computational  examples  of  active  flow  control  simulations  have  been  provided  to  illustrate  utility 
of  the  scheme.  These  include  suppression  of  acoustic  resonance  in  supersonic  cavity  flowfields,  the  control 
of  the  leading-edge  vortex  of  a  delta  wing,  the  enhancement  of  efficiency  for  a  transitional  highly  loaded 
low-pressure  turbine  blade,  and  separation  control  for  a  wall-mounted  hump  model.  Active  control  devices 
considered  in  these  simulations  consisted  of  both  steady  and  pulsed  blowing,  plasma  actuators,  vortex  gen¬ 
erating  jets,  continuous  suction,  and  oscillatory  blowing/suction.  The  variety  of  these  mechanisms  serve  to 
indicate  a  general  applicability  of  the  methodology. 

A  number  of  “pacing”  items  which  require  additional  development,  and  that  can  extend  the  potential 
of  the  scheme  to  other  applications,  are  enumerated  here.  The  most  problematic  of  these  is  the  ability 
of  simulations  to  capture  shock  waves,  and  still  maintain  both  stability  and  high-order  spatial  accuracy. 
Although  the  external  flowfield  of  the  cavity  computation  presented  above  was  supersonic,  the  Mach  number 
was  sufficiently  low  (Mao=1.19)  that  it  did  not  adversely  impact  the  high-order  techniques.  Large-eddy 
simulations  of  supersonic  compression  ramp  flows  carried  out  by  Rizzetta  and  Visbal[39]  with  Mao= 3.0 
employed  a  hybrid  shock  capturing  approach.  Although  this  technique  worked  well  for  a  single  shock  and 
a  two-dimensional  geometry,  it  was  not  deemed  to  be  flexible  enough  for  more  general  applications.  Recent 
investigations  by  Visbal  and  Gaitonde[81]  and  by  Croker  and  Gaitonde[82]  have  explored  alternative  methods. 

The  current  implementation  for  communicating  information  at  overset  boundaries  is  by  a  pre-processing 
operation  to  obtain  interpolation  stencils.  This  is  true  for  situations  where  data  is  exchanged  by  direct 
injection  which  is  used  to  accommodate  domain  decomposition  for  parallel  processing,  as  well  as  the  more 
general  circumstance  of  overlapping  topologically  distinct  grids.  The  ability  to  represent  translating  and/or 
deforming  meshes  that  alter  their  donor  and  receiver  stencils  during  the  course  of  a  simulation,  is  precluded 
by  this  approach.  One  of  the  features  considered  for  future  development  is  the  capability  of  “on  the  fly” 
interpolation,  that  would  allow  simulation  of  relative  motion  between  bodies. 

Finally,  it  may  also  be  possible  to  provide  an  adaptive  component  to  the  numerical  scheme,  and  in 
particular  to  the  spatial  filter.  This  would  permit  the  order  of  the  filter  to  increase  locally,  and  enhanced 
resolution  which  responds  as  features  in  the  solution  evolve.  And  unlike  adaptive  meshing,  no  improved 
grids  would  need  to  be  generated  as  part  of  the  process. 


25 


Figure  1:  Dispersion-error  characteristics  of  various  spatial  discretizations  for  one-dimensional  periodic  func¬ 
tions. 
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Figure  2:  Dissipation-error  characteristics  of  various  filters  for  one-dimensional  periodic  functions. 
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Figure  3:  Dissipation-error  characteristics  with  various  values  of  the  implicit  coefficient  a/  and  a  lOth-order 
filter  for  one-dimensional  periodic  functions. 
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Figure  4:  Schematic  illustration  of  five-point  mesh  overlap. 
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Figure  5:  Schematic  representation  of  plasma  actuator. 
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Figure  6:  Geometry  for  the  empirical  plasma- force  model. 
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Figure  9:  Spanwise-averaged  time-mean  pressure  coefficient  contours. 
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Figure  11:  Spanwise-averaged  time-mean  steamwise  velocity  contours. 


Figure  12:  Spanwise-averaged  time-mean  steamwise  velocity  profiles. 
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Figure  13:  Spanwise-averaged  time-mean  turbulent  kinetic  energy  contours. 
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Figure  14:  Time-mean  turbulent  kinetic  energy  spanwise  wave-number  spectra  at  x  =  0.2, 0.5, 0.8  and 
y  =  0.2  S0/l. 
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Figure  15:  Spanwise-averaged  time-mean  fluctuating  pressure  contours. 
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Figure  16:  Spanwise-averaged  time-mean  fluctuating  pressure  profiles. 
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Figure  17:  Instantaneous  Mach  number  contours  at  the  midspan  location. 
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Figure  18:  Instantaneous  spanwise  vorticity  contours  at  the  midspan  location  for  t  =  t\. 
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Figure  19:  Instantaneous  spanwise  vorticity  contours  at  the  midspan  location  for  t  =  t-2- 
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Figure  20:  Instantaneous  spanwise  vorticity  contours  at  the  midspan  location  for  t  =  £3. 
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Figure  21:  Instantaneous  spanwise  vorticity  contours  at  the  midspan  location  for  t  =  1 4. 
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Figure  22:  Spanwise-averaged  turbulent  kinetic  energy  frequency  spectra  at  x  =  0.2, 0.5, 0.8  and  y  =  0.2Sq/1- 
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Figure  23:  Instantaneous  total  pressure  coefficient  contours  and  iso-surface. 
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Figure  24:  Spanwise-averaged  fluctuating  pressure  frequency  spectra  on  the  cavity  rear  bulkhead  at  y  = 
-0.04. 
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Figure  25:  Spanwise-averaged  fluctuating  pressure  frequency  spectra  on  the  cavity  floor  at  x  =  0.2. 
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Figure  26:  Spanwise-averaged  fluctuating  pressure  frequency  spectra  on  the  cavity  floor  at  x  =  0.5. 
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Figure  27:  Spanwise-averaged  fluctuating  pressure  frequency  spectra  on  the  cavity  floor  at  x  =  0.8. 
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Figure  28:  Schematic  representation  of  delta  wing  actuator  locations  and  flow  features. 
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Figure  29:  Time  history  of  vortex  breakdown  location  for  i?e=25,000  and  AOA=34  deg. 
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Figure  30:  Instantaneous  planar  contours  through  the  vortex  core  for  i?e=25,000  and  AOA=34  deg:  (a)  and 
(d)  pressure  coefficient,  (b)  and  (e)  streamwise  velocity,  (c)  and  (f)  streamwise  vorticity. 
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Figure  31:  Instantaneous  planar  contours  for  i?e=25,000  and  AOA=34  deg:  (a)  streamwise  vorticity,  (b) 
streamwise  velocity. 
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Figure  32:  Instantaneous  iso-surfaces  for  i?e=25,000  and  AOA=34  deg:  (a),  (b),  (d),  and  (e)  streamwise 
vorticity,  (c)  and  (f)  total  pressure. 
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Figure  33:  Time-mean  flow  quantities  for  Re= 25,000  and  AOA=34  deg:  (a)  and  (d)  planar  contours  of 
streamwise  velocity  through  the  vortex  core,  (b)  and  (e)  planar  contours  of  rms  streamwise  velocity  fluctu¬ 
ations  through  the  vortex  core,  (c)  and  (f)  iso-surface  of  streamwise  vorticity. 
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Figure  34:  Instantaneous  planar  contours  through  the  vortex  core  for  i?e=25,000  and  AOA=38  deg:  (a)  and 
(c)  streamwise  velocity,  (b)  and  (d)  streamwise  vorticity. 
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Figure  35:  Instantaneous  iso-surfaces  at  AOA=38  deg:  (a)  and  (c)  streamwise  vorticity  for  i?e=25,000,  (b) 
and  (d)  total  pressure  for  i?e=9,200. 


60 


Figure  36:  Schematic  representation  of  the  turbine  blade  configuration. 
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Figure  37:  Turbine  blade  computational  mesh  system. 
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Figure  38:  Duty  cycle  amplitude  time  history. 
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Figure  39:  Vortex  generator  jet  geometry  and  mesh  system. 
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Figure  40:  Time-mean  surface  pressure  coefficient  distributions  for  vortex-jet  control. 
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Figure  41:  Time-mean  planar  contours  for  vortex-jet  control:  a)  streamwise  velocity,  b)  spanwise  vorticity. 


66 


Figure  42:  Time-mean  turbulent  kinetic  energy  spanwise  wave-number  spectra  for  vortex-jet  control. 
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Figure  43:  Instantaneous  planar  contours  for  vortex-jet  control:  a)  streamwise  velocity,  b)  spanwise  vorticity. 
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Figure  44:  Turbulent  kinetic  energy  frequency  spectra  for  vortex-jet  control. 
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Figure  45:  Instantaneous  iso-surfaces  of  vorticity  magnitude  in  the  trailing-edge  region  for  vortex-jet  control. 
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Figure  47:  Time-mean  surface  pressure  coefficient  distributions  for  plasma  control. 
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Figure  48:  Time-mean  results  for  plasma  control:  contours  of  streamwise  velocity  (top  row),  streamlines 
(middle  row),  contours  of  Cp  (bottom  row). 
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Figure  49:  Time-mean  turbulent  kinetic  energy  spanwise  wave-number  spectra  for  plasma  control. 
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Figure  50:  Instantaneous  streamlines  at  the  midspan  for  plasma  control  with  configuration  A. 
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Figure  51:  Instantaneous  contours  for  plasma  control:  streamwise  velocity  (top  row),  spanwise  vorticity 
(middle  row),  spanwise  vorticity  on  the  blade  surface  (bottom  row). 
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Figure  52:  Turbulent  kinetic  energy  frequency  spectra  for  plasma  control. 
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Figure  53:  Instantaneous  iso-surfaces  of  vorticity  magnitude  colored  by  streamwise  velocity  in  the  trailing- 
edge  region  for  plasma  control. 
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Figure  54:  Computational  mesh  system  for  wall-mounted  hump. 
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Figure  55:  Time-mean  surface  pressure  coefficient  distributions. 
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Figure  56:  Time-mean  streamwise  velocity  contours  for  the  baseline  case:  LES  (top),  RANS  (middle), 
experiment  (bottom). 
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Figure  57:  Time-mean  streamwise  velocity  contours  for  the  suction  case:  LES  (top),  RANS  (middle),  exper¬ 
iment  (bottom). 
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Figure  58:  Time-mean  streamwise  velocity  contours  for  the  oscillating  blowing/suction  case:  LES  (top), 
RANS  (middle),  experiment  (bottom). 
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Figure  59:  Time-mean  Reynolds  stress  contours  for  the  baseline  case:  LES  (top),  RANS  (middle),  experiment 
(bottom). 


84 


0.2 


Figure  60:  Time-mean  Reynolds  stress  contours  for  the  suction  case:  LES  (top),  RANS  (middle),  experiment 
(bottom). 
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Figure  61:  Time-mean  Reynolds  stress  contours  for  the  oscillating  blowing/suction  case:  LES  (top),  RANS 
(middle),  experiment  (bottom). 


86 


Figure  62:  Phase-averaged  spanwise  vorticity  contours  at  0=0  deg  for  the  oscillating  blowing/suction  case 
LES  (top),  RANS  (middle),  experiment  (bottom). 
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Figure  63:  Phase-averaged  spanwise  vorticity  contours  at  0=90  deg  for  the  oscillating  blowing/suction  case: 


LES  (top),  RANS  (middle),  experiment  (bottom). 


88 


Figure  64:  Phase-averaged  spanwise  vorticity  contours  at  0=180  deg  for  the  oscillating  blowing/suction 
case:  LES  (top),  RANS  (middle),  experiment  (bottom). 
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Figure  65:  Phase-averaged  spanwise  vorticity  contours  at  0=270  deg  for  the  oscillating  blowing/suction 
case:  LES  (top),  RANS  (middle),  experiment  (bottom). 
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Figure  66:  Instantaneous  contours  of  streamwise  velocity  and  iso-surface  of  vorticity  magnitude  for  the 
oscillating  blowing/suction  case. 
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=  coefficients  of  explicit  terms  in  the  compact  difference  formula 
=  coefficients  of  explicit  terms  in  the  compact  filter  formula 
=  duty  cycle  amplitude  function 
=  nondimensional  turbine  inter-blade  spacing,  4/3 
=  pressure  coefficient 

=  integrated  wake  total  pressure  loss  coefficient 
=  cavity  depth 
=  plasma  scale  parameter 
=  electron  charge,  1.6  x  10“ 19  coulomb 
=  nondimensional  electric  field  vector 
=  total  specific  energy 

=  nondimensional  turbulent  kinetic  energy  wave  number  and  frequency  spectra 
=  reference  electric  field  magnitude 

=  nondimensional  components  of  the  electric  field  vector 
=  dimensional  imposed  actuator  pulsing  frequency,  Hz 
=  inviscid  vector  fluxes 
=  viscous  vector  fluxes 

=  grid  indices  of  computational  coordinates  £,  77,  / 

=  indices  of  donor  point  stencil 
=  Jacobian  of  the  coordinate  transformation 
=  characteristic  length  for  nondimensionalization 
=  Mach  number 

=  nondimensional  static  pressure 
=  nondimensional  inlet  and  outlet  total  pressure 
=  Prandtl  number,  0.73  for  air 
=  nondimensional  charge  density 
=  turbine  blade  inflow  velocity  magnitude 
=  vector  of  dependent  variables 
=  components  of  the  heat  flux  vector 
=  reference  Reynolds  number,  PooUool/ Moo 
=  Ri ,  Rj ,  Rk  -  interpolation  coefficients 
=  nondimensional  local  span 

=  scaled  spatial  coordinate  for  wave  number  analysis 
=  source  vector 
=  nondimensional  time 

=  portion  of  fundamental  period  over  which  actuator  is  active 
=  nondimensional  actuator  fundamental  period 
=  nondimensional  static  temperature 
=  spectral  transfer  function 

=  nondimensional  Cartesian  velocity  components  in  the  x,y,z  directions 
=  U,  V,  w 

=  contravariant  velocity  components 

=  nondimensional  Cartesian  coordinates  in  the  streamwise,  vertical,  and  spanwise  directions 
=  x,y,z 

=  coefficients  of  implicit  terms  in  the  compact  difference  formula 
=  coefficients  of  implicit  terms  in  the  compact  filter  formula 
=  turbine  inflow  and  outflow  blade  angles 
=  specific  heat  ratio,  1.4  for  air 
=  Kronecker  delta  function 

=  2nd-order  explicit  and  nth-order  compact  finite-difference  operators  in  £,  77,  £ 

=  boundary-layer  thickness 
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=  Qp+1  -  Qp 

=  Aj,  Aj,  A*;  -  interpolation  offsets 
=  time  step  size 
=  spatial  step  size 

=  boundary-layer  momentum  thickness 
=  angle  for  phase  averaging 
=  nondimensional  molecular  viscosity  coefficient 
=  nondimensional  body-fitted  computational  coordinates 
=  metric  coefficients  of  the  coordinate  transformation 


=  nondimensional  fluid  density 
=  electron  charge  number  density,  1  x  1011/cm3 
=  oder-of-accuracy  of  interpolation  formula 
=  components  of  the  viscous  stress  tensor 

=  general  function 

=  dimensional  or  nondimensional  frequency 
=  scaled  spatial  wave  number 
=  modified  scaled  spatial  wave  number 


Subscripts 

b 

i 
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max 

min 

P 
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oo 

Superscripts 

n 

P 


=  breakdown  location 
=  grid  point  index  number 
=  Fourier  mode  number 
=  maximum  value 
=  minimum  value 
=  receiving  interpolation  point 
=  numerical  approximation 
=  dimensional  reference  value 

=  time  level 
=  subiteration  level 
=  Fourier  component 
=  d/dS 

=  filtered  value 
=  interpolated  value 
=  time-mean  quantity 
=  fluctuating  component 
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